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TECHNICAL MEMORANDUM 


HIGH PERFORMANCE, ROBUST CONTROL OF FLEXIBLE SPACE STRUCTURES 
(MSFC Center Director’s Discretionary Fund Final Report, Project No. 96—23 ) 


1. INTRODUCTION 


Where there is no vision, the people perish. Proverbs 29:18. 

Throughout history, mankind’s ambition has fueled a drive to explore. The last four decades 
in particular have seen a great expansion of technology brought about by the space program. Along 
with advances in technology comes the opportunity to expand our understanding of the universe which 
demands even greater advances in technology. To realize ambitious technical goals requires increasingly 
sophisticated spacecraft systems. One technology area that is especially challenging is pointing control 
systems for flexible space structures. Applications include space telescopes and interferometers. Earth- 
observing spacecraft, communications spacecraft, and military spacecraft such as space-based lasers 
and tracking systems. Control systems for flexible space structures must meet stringent performance 
requirements while being robust to model uncertainties. The topic of this dissertation is robust control 
system design for flexible space structures (FSS’s) that achieve high performance on orbit. 

1.1 Challenges of Flexible Space Structure Control 

Vibration control of flexible spacecraft has motivated intense research in modern and robust 
control theory. In general, FSS's possess several pathologies that compound the difficulty of controller 
design. Light-weight, low-stiffness structures coupled with stringent pointing requirements result in 
numerous closely spaced, lightly damped vibrational modes within the controller bandwidth. A primary 
difficulty with control of FSS’s is that accurate models of the on-orbit system are not typically available 
a priori. Successful implementation of a high-performance, on-orbit pointing control system requires 
integration of system identification and control system design in a form that facilitates iterative redesign, 
or “tuning,” of the control system. 

Recent experience with the Hubble Space Telescope (HST) serves to illustrate the need for on- 
orbit system identification and controller tuning. 2 After initial deployment, the HST experienced severe 
pointing disturbances due to thermal excitation of the solar arrays and the performance degradation was 
attributed to model errors. On-orbit test data were used to update the control design model for controller 
redesign. Incremental control augmentation in the form of forward path and inner loop feedback filters 
was employed with each subsequent retune phase exhibiting increasing levels of performance. However, 
this redesign effort took considerable time and a very substantial investment of manpower. Considerable 
savings could have been gained by a more efficient methodology for on-orbit control system redesign. 



For many spacecraft, to obtain high-performance pointing control on orbit requires closed-loop 
system identification and controller refinement to improve performance. A methodical process for 
on-orbit tuning for flexible spacecraft must consist of Five steps: ( I ) Empirical or analytical model 
development for initial control design; (2) low authority control design for robust stability; (3) closed- 
loop system identification to improve knowledge of the on-orbit plant dynamics; (4) controller refine- 
ment to increase control authority and improve performance; and (5) performance assessment. 

Steps (3 ) — (3) may be iterated until the desired controller performance is obtained. 

Closed-loop system identification (step 3) is significant for a number of reasons. Typically, 
spacecraft cannot be tested on orbit in the open-loop manner. Also, open-loop modeling errors are not 
necessarily indicative of closed-loop robustness or performance. 3 For control system tuning with flexible 
spacecraft, the plant model should be obtained from closed-loop excitation in the orbital environment. 

A new procedure for closed-loop system identification is also presented which tunes the parameters 
of an open-loop (state space) control design model to better characterize the on-orbit dynamics. 

Recent advancements in the field of robust multivariable control have provided a fruitful para- 
digm for on-orbit controller (step 4) tuning where the competing issues of performance and stability 
are explicitly addressed. Mixed H 2 I H w control theory is well suited to on-orbit controller tuning since 
the technique generates a set of controllers that trade between robustness and performance. The control 
formulation is such that the controller parameters are tuned to account for plant changes or varying 
performance requirements. Successive implementation of these controllers allows the compensator 
to be selected which achieves maximum performance in the presence of model errors. In this research 
effort, a homotopy algorithm was developed for synthesis of fixed-order, mixed H 2 tH ao compensators. 

1.2 Contributions 

It is the goal of this dissertation to develop a control design and system identification procedure 
that achieves high performance while providing robust stability guarantees. Contributions of this 
research effort include the following: 


1. This dissertation will present a homotopy algorithm for fixed order, mixed H 2 / 
compensator synthesis. 

2. To successfully implement this homotopy algorithm required the development of a new 
algorithm for numerical optimization that is robust to ill-conditioned and indefinite Hessians. The utility 
of this new optimization algorithm is demonstrated by examples in appendix A. 

3. Examples are given that illustrate the significance of mixed-norm design for nominal perfor- 
mance and robust stability. Included is a control design exercise for an experimental flexible space 
structure that includes open-loop system identification for control design model development. 

4. A new method for system identification from closed-loop response data is presented for 
multi-input, multi-output systems in canonical form. Examples are given to demonstrate the closed-loop 
system identification and controller redesign process. 
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This dissertation is structured as follows. In section 2, a brief survey of the developments in 
control design and system identification which relate to controller tuning will be presented. Special 
attention is given to fixed-order control design for both robustness and performance. Emphasis is also 
placed on recent developments in closed-loop system identification. The approaches to control design 
for nominal performance, robust performance, and nominal performance/robust stability are highlighted 
in section 3. The formulation of the fixed-order mixed H 2 I H ^ control design method is also given. 
Homotopy methods are introduced in section 4 and a homotopy algorithm developed to synthesize 
mixed H 2 ///«, controllers is presented. Section 5 presents mixed-norm control design examples. 

The new closed-loop system identification procedure is given in section 6 with examples. Conclusions 
and recommendations are provided in section 7 and appendix A is included to present the details of the 
numerical optimization procedure developed for ill-conditioned and indefinite Hessians. 
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2. BACKGROUND 


2.1 Fixed-Order Control Design For Stability and Performance 

Recent years have seen significant advances in modem control theory with H 2 and H x methods 
gaining widespread recognition and application. Early work in multivariable control synthesis utilized a 
quadratic cost functional to minimize the H 2 norm of a system response to white noise inputs. While the 
H 2 procedure is well suited to design for nominal performance in terms of root-mean-square (RMS) 
quantities such as minimizing line-of-sight errors or control energy, it is well known that stability and 
performance cannot be guaranteed in the presence of model uncertainties. Robustness is addressed in 

control theory which guarantees stability and performance (when defined by an H ^ norm measure) 
in the presence of unstructured uncertainty models, albeit often resulting in overly conservative designs. 
H -analysis and jj -synthesis methods have reduced the conservatism of methods by accounting 
for the structure in the uncertainty and performance specifications. 4-6 

A significant disadvantage of these modern control techniques is that the resulting compensator 
is the same order as the generalized plant, which is often larger than the original plant due to the inclu- 
sion of frequency-dependent weights to achieve the desired performance and robustness characteristics. 
The consequential large controller order can be indirectly alleviated by reducing the order of the control- 
ler or alternatively by reducing the order of the design plant. In either case, indirect methods are subopti- 
mal in performance and do not guarantee closed-loop stability. However, direct methods may be em- 
ployed which impose constraints on controller order or architecture in the optimization procedure 
and hence provide stability and performance guarantees. In an optimization-based synthesis procedure, 
necessary conditions are formulated for the constrained closed-loop system that ensure internal stability. 

The Optimal Projection method of Hyland and Bernstein is an H 2 procedure whereby order 
constraints are imposed on the controller and the necessary conditions for minimizing a quadratic cost 
functional with respect to the fixed-order controller are derived. 7 The resulting necessary conditions 
consist of two modified Riccati equations and two modified Lyapunov equations coupled by an oblique 
projection matrix. However, solution of the necessary conditions for realistic large-order systems is a 
difficult task. Homotopy methods have been employed to solve the optimal projection equations. s 

As a means of providing robust stability with RMS-type performance, the mixed // 2 / 
methodology has been developed. Much work has been done with variations of the mixed H-, / 
problem (for a summary see reference 9). The earliest refereed work was done by Bernstein and Haddad 
who extended the Optimal Projection approach to (linear-quadratic-Gaussian) LQG control with an H 2 
norm constraint. 1 Their formulation minimized an overbound on the H 2 norm from a disturbance input 
to one output while satisfying an H x norm overbound from the same disturbance input to a second 
output. This approach is potentially conservative as a result of minimizing an overbound on the H 2 
norm. This approach was extended to the general mixed H 2 / H x problem which can be specialized to 
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both full- and fixed-order pure H 2 and pure H x controllers.*^ Reference 1 1 further extended the mixed 
Hi / //oo formulation of reference 1 by fully accounting for singularities in the problem formulation. 
Rotea and Khargonekar formulated the general mixed Hi ! H x problem with independent inputs and 
outputs for the two transfer functions and minimized the actual Hi norm based on full state feedback. 12 
Ridgely extended the formulation to output feedback including the fixed-order case with either regular 
or singular H x constraints 9 and provided a numerical solution.* 3 Another approach to the general 
mixed Hi / H^ problem was developed by Sweriduk and Calise 14 who used a differential game formu- 
lation with a conjugate gradient algorithm to obtain fixed-order controllers. In this dissertation, a nu- 
merical solution of this formulation using a homotopy algorithm is developed and presented in section 4. 
Based on the formulation of reference 1 , homotopy algorithms have also been recently developed for 
both discrete time and continuous time mixed H-> / H x design.* 5 -* 6 

Ridgely et al. shows the mixed Hi / H problem to be a strictly convex optimization problem 
with a unique solution when the controller is of the order equal to or larger than the underlying // 2 
problem. 9 The solution is shown to lie on the boundary of the infinity norm constraint when active 
(for y < y where y is the H ^ norm when the optimal Hi controller is used) and is just the Hi con- 
troller when y > y . For controllers with order less than the underlying Hi problem, the possibility of 
local minima in the unconstrained Hi problem (Optimal Projection) exists and the solution of the mixed 
Hi / H^ problem may or may not lie on the boundary of the H ^ norm constraint. As will be shown 
in the next section, Sweriduk and Calise 14 begin with the fixed-order H^ cost functional and append 
the fixed-order // 2 cost functional. Reference 9 takes the converse approach by appending the H ^ norm 
constraint to the // 2 cost functional. As a consequence, the formulation in reference 9 can handle both 
regular and singular H^ constraints and can be specialized to the Hi problem. Conversely, the refer- 
ence 14 formulation can handle both regular and singular Hj constraints and can be specialized to the 
H problem. However, the ability to handle singular constraints on either cost exists as long as one 
of the cost constraints is regular. Whereas the reference 14 formulation assumes that the plant dynamics 
(A matrix) of the two transfer functions are the same, reference 9 allows different dynamics in the two 
separate underlying H ^ and Hi problems. Using the canonical compensator with a static gain output 
formulation presented in section 3, the reference 14 formulation requires the simultaneous solution 
of five coupled nonlinear matrix equations. Reference 9 presents the necessary conditions in the form 
of seven coupled, nonlinear matrix equations. 

2.2 Closed-Loop System Identification and Control Design 

Implicit in the design of high-performance control systems is the availability of an accurate 
model of the system to be controlled. Although system identification and control design are both critical 
aspects of high-performance model based control design, the theoretical foundations of these two disci- 
plines have developed distinctly. Developments in system identification have been directed toward 
obtaining accurate nominal models with bounds on the associated uncertainty. Recognizing the depen- 
dence on accurate nominal design models, robust control theory has been developed to accommodate 
modeling errors. Recently, attention has been drawn to the fact that the issues of system identification 
and control design must be treated as mutually dependent. In reference 3, Skelton points out that since 
the magnitude and spectrum of excitation forces are controller-dependent, an appropriate model 
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for control design cannot be determined independent of the controller. The point is made that the validity 
of the model is dictated by the controller instead of the opposite, as is usually assumed. Since the model 
that is most appropriate depends on the control design, the open-loop response of a model is not suffi- 
cient to indicate the fidelity of the model for control design. Additionally, robust performance requires 
an accurate model of the plant in the controller crossover frequency range, 17 indicating that the amount 
of model error that can be tolerated is frequency- and controller-dependent. Hence, the issues of model 
identification and model-based control design must be treated as a joint problem suggesting an iterative 
solution.- 1 • lx 

Closed-loop system identification (i.e., identification of the open-loop plant given closed-loop 
response data and knowledge of the compensator dynamics) is currently a field of active research. 

Most of the methods that identify state space models of the open-loop plant are based on identifying 
the closed-loop Markov parameters from frequency response functions and then extracting the open-loop 
Markov parameters. A discrete time state space realization is then obtained from these identified Markov 
parameters. An approach to iterative closed-loop system identification and controller redesign is given 
by Liu and Skelton where a state space model is estimated using the q-Markov Cover algorithm and 
a controller is designed using the Output Variance Constraint algorithm. 19 q-Markov Cover theory 
describes all realizations of a linear system which match the first q-Markov parameters and covariance 
parameters of the true system generated by pulse responses. First, the entire closed-loop system is 
estimated and then, using knowledge of the compensator dynamics, the plant is extracted. However, the 
identified open-loop plant has dimension equal to the controller dimension plus the dimension of the 
identified closed-loop plant. Model order reduction is used to remove the superfluous states. A similar 
approach is taken in reference 20 where Phan et al. formulate a method which first obtains the closed- 
loop Markov parameters and then the open-loop Markov parameters are recovered. A discrete time state 
space realization is then obtained from the Markov parameters. Other related approaches are given 
in references 2 1 and 22. 

To provide synergism, an iterative process should match the system identification objectives 
with the control design objectives. Reference 18 presents an iterative algorithm for frequency-response 
identification from closed-loop data and robust control design. The identification phase is control- 
oriented with the objective of providing robust performance by closely approximating the achieved 
closed-loop performance. The interplay between identification and control design is formalized 
by specifying a performance metric (norm) for model-based optimal control design and an identification 
cost function that minimizes the difference in the achieved performance and nominal design perfor- 
mance as defined by the same metric. Many algorithms have been developed based on this general 
framework which utilize different performance norms and identification algorithms. An excellent survey 
of this topic is given in reference 23. Most of the work cited therein has been done in the context 
of single-input, single-output (SISO) systems. 

A classical approach to parameter estimation which has been extended to closed-loop system 
identification is the prediction error approach. 24 This optimization method estimates the parameters 
of a linear system by minimizing the squared sum of the errors between the actual measurements 
and the predicted measurements. Zang, Bitmead, and Gevers present an iterative prediction error 
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identification and control design algorithm based on the Hj norm. 2 - s The control objective is used 
to frequency-weight the identification cost functional and the resulting prediction error spectrum is used 
to frequency-weight the control design. In this manner the control is penalized heavily where the SISO 
transfer function model fit is poor and the model is weighted to fit best in the regions most critical 
to performance. Another approach utilizes the dual Youla parameterization of all plants stabilized 
by a given controller. This approach was introduced by Hansen and Franklin-^ and further elaborated 
by Hansen et al. in reference 27 as applied to closed-loop experiment design. Schrama applied the dual 
Youla parameterization to closed-loop system identification in references 28 and 29 as did Anderson 
and Kosut in reference 30. A related method directed toward identification of the coprime factors 
of the plant was introduced by Schrama in reference 28 and was further elaborated in references 29 
and 3 1 . 


In this dissertation the prediction error method is extended to closed-loop identification for 
multivariable systems. Building on a method developed for estimation ol the parameters of an open-loop 
system in canonical form from open-loop data, a new procedure for closed-loop system identification 
is developed and demonstrated in section 6. 
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3 . APPROACHES TO CONTROL DESIGN FOR FLEXIBLE SPACE STRUCTURES 


H 2 methods are often used when designing control systems to reduce the vibration response 
of a flexible structure. While H 2 design gives good nominal performance, the controllers are highly 
tuned to the design model and errors in the design model are not accounted for, typically inducing 
instability at higher levels of control authority. As a result, the actual performance achievable is limited 
with H 2 designs. To achieve high levels of performance in the actual system, robustness to model errors 
must be taken into account in the design process as is done in the H M control design theory. In the 
following sections, a brief introduction to H 2 ,H 0 Q . ju -synthesis, and mixed H 2 //j control theory 
is given followed by the problem formulation for fixed-order control design. 

3.1 Control Design for Stability and Performance 

The generalized plant of a standard control problem is given by 


.V = Av + By w 4- 

in 

Z = Cy X + Dy jli 

(2) 

— C"2 A* + D~> 1 vv + Z? 22 ^ 1 

(3) 


where x e R" is the state vector, vv s R' m is the disturbance vector, u e R nu is the control vector, 
r e R"' is the performance vector, and y € R ny is the measurement vector. The following is assumed: 

• (A, 5 |, Cj) is stabilizable and detectable 

• D|2 has full column rank 

• D 2 1 has full row rank. 

A general compensator for this system is 


-V c ®f.v (4) 

a = C c x c , (5) 

where „v t . e R n< is the state vector of the controller, the dimension of which can be specified. Closing 
the loop using negative feedback yields the closed-loop system dynamics 


-V = Ay + Bw 


( 6 ) 
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where 



( 8 ) 


A - B~,C r 

A = 

B^C'i A ( . — B c Dr2C c 

- r B \ 1 

B = 

_ ^( ^21 

c = (c, -z> l2 c ( . ] 

The set of all internally stabilizing compensators is defined as 


(9) 


( 10 ) 

(ID 


S c = j(/4 r , B c ,C c ) :A is asymptotically stablej . (12) 

3.1.1 Nominal Performance Design 

For an Hj problem, the objective is to minimize the Hj norm of the closed-loop transfer 
function from disturbance inputs to performance outputs 


T = 

1 rvv 




(13) 


where the disturbances are confined to the set of signals with bounded power and fixed spectra. If the 
disturbance is modeled as white noise, the objective is 


min 



lim E 

/— »oo 


{--(<) r --(<)} 



(14) 


where £{*} is the expectation operator. The cost can be expressed as 

j(A c ,B c ,C c ) = trj<25B r J = trjpc^cj , (15) 

where 

AP+PA 7 +BB T =0 (16) 

A T Q + QA+C T C = 0 . 


(17) 



P is the controllability grammian of (-4.5) and Q is the observability grammian of (c./i). An equiva- 
lent cost functional also arises for the case of impulsive inputs. 

Another approach to design for nominal performance employs the H x norm, which can be 
interpreted as the gain of the system and is the worst-case amplification over all inputs vc (/) of bounded 
energy Ln signals. From a frequency domain perspective, the H ^ norm is defined as the maximum 
singular value of T(s) over all frequencies, i.e.. 

ll 7 'rv, IL=^p{a(r rM .(./ro))} ^ (18) 

CO 

where o denotes the maximum singular value. H x control design theory, based on references 32 
and 33, involves defining (possibly frequency-dependent) weights on the inputs and outputs such that 
the performance objectives are satisfied by minimizing ||7’; M >|| oo - Because the norm is defined with 
respect to the peak magnitude of the transfer matrix frequency response and the H 2 norm is defined 
by an integral square quantity (in time or frequency by Parseval’s Theorem), the respective closed-loop 
systems may have considerably different characteristics. Depending on the performance objectives, one 
design procedure may be preferable to the other. When performance is specified by RMS measures, H~, 
design typically yields better nominal performance. The significant benefit of H ^ theory is that robust- 
ness to unstructured model errors is explicitly factored into the design process. 

3.1.2 Robust Stability and Robust Performance 

In addition to nominal performance, robust stability is an important design consideration. Robust 
stability requires the closed-loop system to remain stable for bounded model errors. The uncertainty may 
be modeled in many forms such as multiplicative, inverse multiplicative, additive, parametric, etc., and 
may be located at various points in the loop. By absorbing all of the scalings and weights into the plant, 

P , the robust stability problem may be formulated as the linear fractional transformation (LFT) in 
figure 1 . The uncertainties are scaled so that A^ is the set of all stable perturbations, A , such that 
HAL < 5 . Assuming that K(s ) internally stabilizes the closed-loop for A = 0, then a sufficient condition 
for robust stability for all plants in the set formed by A e A^ is that 34 - 35 

• ( i9) 

Thus, like the nominal performance problem, robust stability is characterized by the H norm 
of a transfer function. 

It is the ability to formulate the performance problem as a robust stability problem that enables 
design for robust performance in the H x setting. Consider the LFT of an uncertain plant in figure 2 
with inputs and outputs defined for performance and an uncertainty model. The conditions for robust 
performance are: 

1. Robust stability (equation (19)) 

2. Performance maintained for all AeA^. 
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w 
u 

Figure I . LFT for robust stability analysis. 


A 



Figure 2. LFT for robust performance design. 



Closing the loop from r 2 to w 2 through a fictitious uncertainty block A p recasts the robust performance 
problem as a robust stability problem, shown in figure 2, where the blocks are scaled to one. 

A sufficient condition for robust performance is that 

l|7WL<|. (20) 

Define to be the set of all stable, bounded, unstructured perturbations A such that flAL < S. When 
A e A^, equation (20) is necessary and sufficient to ensure robust stability. Designing for robust perfor- 
mance using A p as in figure 2 introduces a block diagonal structure to A which results in equation (20) 
being only sufficient and possibly overly conservative. This conservatism is relaxed in the /J -analysis 
and fj -synthesis procedures by accounting for the block diagonal structure in A. The structured singu- 
lar value, fi , is the inverse of the smallest destabilizing perturbation of a transfer function matrix and is 
defined as: 


v(T-.w{jM)) = 


1 

m//;|d r (A(y'tw)) : A e A^detf/ - T :w (j(o)A(jco)] = 0} 


( 21 ) 
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The structured singular value is used to define the /u -measure, which although not a norm, is given by 


P’U (0 t ti = *up{K T :*-U a >))} • ( 22 ) 

Hence, the sufficient condition for robust performance in equation (20) becomes the necessary 
and sufficient condition 

P'( ja> % < ^ ■ (23) 

Although the structured singular value cannot be directly computed, an upper bound can be computed as 

A/(T)=inf{a(D7Z)- 1 )} . (24 ) 

where D = diagjdy/y j has the same structure as A and dj are scalar, positive, real functions of fre- 
quency. An iterative scheme is used to solve this optimization problem. In the first step, an controller 
is designed and in the second step, the D-scales are optimized for this controller in accordance with 
equation (24). A curve-fitting step is used to determine stable, minimum phase rational transfer function 
fits (in magnitude) to the optimum D-scales. In the next iteration, these D-scales are incorporated into 
the generalized plant and the control design is repeated, followed again by D-scaling. This iterative 
process continues until the upper bound in equation (24) cannot be reduced significantly. 

3.1.3 Robust Stability and Nominal Performance Design 

Although (A -synthesis provides stability and performance in the presence of model errors, the 
performance is defined by an H x norm measure. The mixed H 2 / design procedure has been devel- 
oped to provide robust stability and nominal ( H 2 ) performance by minimizing the H 2 norm for one set 
of inputs/outputs while satisfying an H ^ norm overbound for another set of inputs/outputs. With refer- 
ence to figure 2, the objective is to satisfy 


min 

K 


Z~i w-> 


I? 


(25) 


subject to 


T. 


VTj 


<r . 


( 26 ) 


This problem has been solved for controllers of fixed dimension in references 9, 14, and 36-38. 


The formulation of reference 14 is presented in the next section and a numerical homotopy 
algorithm for this formulation of the mixed H 2 / control design problem is presented in section 4. 
The homotopy algorithm is a two-parameter iterative scheme which effectively trades between robust 
stability and nominal performance by varying the overbound on the H x norm, y , and the weight 
on the H 2 cost, A. For a given y , X is increased until the H ^ norm constraint becomes an active 
equality constraint (at which point the H 2 norm can no longer be reduced), or until the H 2 norm ceases 
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to decrease. The set of controllers where the H x norm is equal to the overbound are called the boundary 
solutions. This set provides an explicit tradeoff between nominal performance and robust stability. By 
incorporating the D-scales from fu -synthesis into the subproblem, the structure of the uncertainty 
block may be accounted for. resulting in a fixed-order mixed H2 I H design procedure. Recent develop- 
ments in robust control theory that avoid this multiplier-controller iteration process include the use of 
absolute stability theory such as Popov analysis . 37-40 This approach provides a means for simultaneous 
optimization of a fixed-architecture controller and a fixed-structure multiplier that accounts for both 
complex and real-structured uncertainty. A homotopy algorithm was developed for synthesis of robust 
controllers with fixed-structure multipliers . 41 

The conflicting demands of designing for robustness and performance may be explicitly ad- 
dressed in the mixed Hj / design setting. The utility of mixed-norm design is exploited by separat- 
ing performance and robustness using the subproblem to design for nominal performance and using 
the H oo subproblem to achieve robust stability. In figure 2 the inputs and outputs associated with the 
uncertainty model comprise vt»| and Z \ , w’2 is associated with disturbance excitations and measurement 
noise, and r 2 is associated with performance outputs. In this manner, the set of boundary controllers 
explicitly trade off nominal performance with levels of robustness to uncertainty. As illustrated in 
figure 3 , this set of controllers includes, at one extreme, the lowest authority controller which provides 
maximum robustness and minimum performance and, at the other extreme, a maximum performance 
controller with the minimum stability margins. Hence, this paradigm is well suited for tuning a control- 
ler on orbit by incrementally increasing control authority. 



Figure 3 . Performance versus robustness trade. 


Another significant aspect of this approach is that each controller in the set may be successively 
implemented to determine the point at which achieved performance begins to deviate from nominal 
performance. At that point, the maximum achievable performance with robust stability for a given 
design model is achieved. The same implementation strategy could be applied to H 2 and (j. control 
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design. However, mixed-norm design tends to allow higher performance levels to be attained for a given 
bounded model error by designing for robust stability. Since H 2 design does not account for model 
errors, the sensitivity to model errors typically results in achieving lower levels of performance than 
the mixed-norm designs. On the other hand, ju designs are not optimized for H 2 performance and 
although robust stability is given, the performance is not comparable to H 2 nominal performance 
design. These issues are illustrated by design examples in section 5. 

In the following section, the problem formulation and necessary conditions for fixed-order 
Hi, H x , and mixed Hi/ H^ control design are presented. Section 4 presents the homotopy algorithm 
for fixed-order, mixed-norm synthesis. 

3.2 Fixed-Order Control Design Formulation 

In order to obtain the H 2 -optimal compensator, the Lagrangian is formed by augmenting the cost 
functional, equation ( 15), with the constraint, equation (17), yielding 

L{Q.LA c ,B c ,C c ) = tr{QBB T +(A T Q + QA + C T Cy} , (27) 

where L is a symmetric matrix of Lagrange multipliers. Matrix gradients are taken to determine the first 
order necessary conditions: 


c]L_ 

dQ 


= 0 . 



dL 

dA c 


= 0 , 



dL 

dC c 


= 0 


(28) 


Hence, computing an # 2 -optimal controller of fixed-order n c < n for the general controller structure 
given in equations (4)-(5) requires the simultaneous solution of five coupled equations. This is not 
only computationally expensive, but is also further complicated by the fact that the problem is over- 
parametrized with such a compensator. 

To optimize the controller for a specified controller state dimension, a “controller form” architec- 
ture is imposed on the compensator dynamics. 42 This minimal realization avoids the problem of over 
parametrization and is a canonical form under mild conditions. 43 The internal structure of the compensa- 
tor is prespecified by assigning a set of feedback invariant indices, v,- . In controller canonical form, 
the compensator is defined as 


\ + n\ - N°y 

(29) 

u c = -Px c 

(30) 

it = — H.x c , 

(31) 
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where v ( . e R IU and n c e R ll} . P and H are free-parameter matrices, and P {) and N () are fixed matrices 
of zeros and ones determined by the choice of controllability indices. Similarly, a compensator in ob- 
server canonical form can be constructed. In this paper only the controller canonical form is employed, 
which imposes the lower bound n c . > ny on the order of the compensator. 


By augmenting the generalized plant of equations ( I )-(3) with the compensator of equations 
(29)-(3 1 ), the system may be written as 



A 

0 " 




b 2 

0 

.V = 

-n°c 2 

p<> 

,v + 

-N {) D 2 1 

VI' H- 

-N°D 2 2 

N° 


= Ax + 5|Ve + B 2 u , (32) 


z — [ Cj 0 j .v + [ D[ 2 0 j u — Cj.v 4 - D\ 2 lt 


(33) 


y = [0 /] .v = C 2 -V , 


(34) 



v = —G v , 


(35) 


where the augmented control input is 


ti = 


u 

u,. 


(36) 


Equations (32)-(35) define a static gain output feedback problem where the compensator is represented 
by a minimal number of free parameters in the design matrix, G. The closed-loop system is given by 


.v = (A — B2GC2 ).v + B\W = Ai + Bw (37) 

r = (c, -D l2 GC 2 )\ = Cx ■ (38) 

Minimizing the H 2 norm of T- w = c|s/ - A j B utilizes the same Lagrangian as given in equation (27), 
but now L is only a function of three parameter matrices, i.e„ L(Q,L,G). Thus, the resulting first order 
necessary conditions are 



( 40 ) 


— = a t q+qa+c t c = o 

dL 

■p 1 = 2(Z)|2 Dj2 GCt - D|2 C] - fij Q)lcT = 0. (41) 

Controller canonical forms can also be used to solve the H ' x problem. The objective is now 
to minimize the norm of the transfer function from disturbance inputs w to performance outputs r 
given by 


T. u 


sup 

\veL-> 



sup ||7: M .|, 

MM' 


(42) 


From reference 14. the optimization problem is to find 

^/ r {||7'ru (G)|| oo : C7 e <3} = y (43) 

where G= | G e )R(nii+ny)xm j s stablej . A more practical design objective is to find a G that ensures 

ll^rwLo - Y f° r some y > y*. This optimization problem may be formulated as a min-max problem 
using the cost functional 


y r KG) = £{/;(/.- ,.)*} , ,44, 

where the expectation operator, £{*}, is taken over a distribution of initial conditions with zero mean 
and variance BBT. In the min-max problem, the minimizing player acts first, so the loop is closed 
around G and the disturbance vv* which maximizes Jy is determined. Then the control G which 
minimizes Jy(w*,G) is determined, resulting in the upper value 


Jy = min max |./ j,(m , ,G)| 
GeGy vveZ^ L J * 


(45) 


G is restricted to the set of admissible controllers defined by Gy = \G e SR(™+m )xm : 

A is stable and ||7'- M ,|| oo < y\. 

For any G e Gy, a unique worst-case disturbance exists and is a feedback of states given by 

n- - y~ 2 B T Q^ , (46) 


16 



where Q x is the positive semidefinite solution of the Riccati equation 

A^ Qoo + QooA + C T C + y ~Q ao BB r Q ao = 0 . (47) 

Substituting the worst-case disturbance from equation (46) into equation (44), the performance index 
becomes 

J y (G) = E^x T (c T C-y- 2 Q 00 BB T Q 00 )m j . (48) 

The objective of minimizing |7'- M ,|| using a fixed-order controller can be fomiulated 


min 

GeGy 


{,/ 7 (G) = tr{(2^5 r }} ^ 


(49) 


subject to equation (47). Using the Lagrange multiplier matrix L, the Lagrangian is 

L{Q 00 ,L,G) = \r\Q 00 BB T + [a T Q 00 + Q x A + C T C + Y~ 2 Q„BB T Q 00 )l} . (50) 

Matrix gradients are taken to determine the first order necessary conditions for an H M suboptimal fixed- 
order compensator gain G: 

-^- = (a + Y~ 2 BB T qJ)L + l(a + Y~ 2 BB T qJ) + BB t = 0 ( 5 1 ) 

= A r Q 00 + Q X A + C T C + Y~ 2 Q°°BB T Q°° = 0 (52) 

oL 

~ = 2( Z)p DpGCi - D^C X - Bl qJlC 7 = 0 . (53) 

oG ' " ' 

As in the Hi problem, three coupled equations have to be solved to obtain a fixed-order compensator 
which satisfies the constraint ||7'_ M .|| < v. 

This formulation minimizes an upper bound on the Hi cost which corresponds to an entropy 
controller. 44 Although this controller guarantees an H ^ norm constraint, as written it is not a pure H 
controller in the sense of reference 33, but can be specialized to the pure H controller. 10 
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Fixed-order H x design has also been extended to fixed-order fj. -synthesis. 36-38 ' 45-47 j^e 
approach taken in reference 36 differs from that of references 45-47 by using optimized fixed-structure 
multipliers instead of multiplier-controller iteration and is based on the earlier developments in refer- 


ences 37 and 38. Because H x controller design is a subproblem when designing for robust performance 


with structured uncertainty, the fixed-order technique just introduced has the potential to constrain the 
order of the controller which is normally subject to significant increases in the jj. -synthesis procedure. 

The mixed // 2 / problem can be approached in a similar fashion. 14 In this case, the general- 

ized plant has additional inputs and outputs with dynamics given by 

a — Ay + BpW'p + Z?|M’ 4- Bin 

(54) 

~P ~ C p ' + D\ pH 

(55) 

Z = C\X + /)[ *> 14 

(56) 

v — C 2 v + & 2 pWp + D~) | vr + D 21 K * 

(57) 

where w p e R»"P and z /} e R n ~ p are the inputs and outputs defining the W 2 subproblem, and w e R uw 
and r e R"~ are the inputs and outputs defining the H ^ subproblem. It is assumed that (4 ,j 3 2 ,C 2 ) 
is stabilizable and detectable. Using the controller canonical form for the compensator, the augmented 
system for the mixed problem is 

A — A A" + Z?| VI ’ + Bp M ' p + Bjll 

(58) 

-p — C p.\ + D\pU 

(59) 

: — C].v + £)| 2 m 

(60) 

y = C 2 .v 

(61) 

u = -Gy , 

(62) 

where 


~ b n 

B P = -N {) D-, p ’ C P = l C P °]’ D ip = [ D \i 

, 0] . (63) 
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The other expressions are the same as in equations (32)-(35). Consequently, the closed-loop system 
is given by 

.v = (A — B 2 GC 2 )■' + BpWp + fi|U' 

= Av - 1 - BpWp + Bw ( 64 ) 


z p — ( C p Di p GC 2 )-v — C p .\ 

(65) 

r = (C, -D i2 GC 2 ) v = C-v . 

(66) 


In order to formulate the performance index of the mixed problem, the Lagrangian for the H 2 problem 
is adjoined to the Lagrangian for the H x problem in equation (50) by a scalar weight, A : 

L = tr [q^BB 7 + AXCj,C p + (A T + Q^A + C T C 

+Y~ 2 QooBB t q x )l + (ax + xa t + B p B T p )l p } . (67 ) 

Note that the H 2 portion of the mixed-norm Lagrangian uses the dual form of the cost and constraint 
in equations ( 1 5) and ( 1 6). This problem was first formulated in reference 1 0. The weight. A . on the H 2 
norm allows a tradeoff between performance ( H 2 norm) and robustness ( norm). The first order 
necessary conditions are 


dL 

dQcc 


T 

= (A + y~~BB T qJ^L + l[a + y~ 2 BB T Qoo) + BB T = 0 


( 68 ) 


^ = A T Q^ + Q^A + C T C + Y^QooBB 7 ^ = 0 
oL 


(69) 


^ = A T Lp+L p A + lC T p C p =0 


(70) 


— = TX + XA T + B b t = 0 
dL p 


(71) 


dL_ 

dG 


[d? 2 D l 2 gc 2 - dJ 2 C, - bJ qJ)L 
+(AdJ D] p GC 2 - XD\ p C p -bJl p )x 1 cj = 0 


(72) 


In section 4, a homotopy method is presented to synthesize fixed-order H 2 , H rjo , and mixed 
H 2 / controllers. 
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4. HOMOTOPY ALGORITHM FOR MIXED H 2 /H 00 CONTROL DESIGN 


Using the formulation and necessary conditions of the fixed-order mixed H 2 ///<*, control design 
problem in section 3. this section presents a homotopy algorithm for synthesis of fixed-order mixed 
H 2 /H oa controllers. Fixed-order H 2 and H x controllers are obtained in the same manner, but for the 
sake of brevity only the mixed-norm algorithm is presented here. A complete development of the H 2 
and H x algorithms is given in reference 48. The approach and algorithms developed here are patterned 
after references 49 and 50. Before embarking on the homotopy algorithm for mixed-norm design, 
an introduction to homotopy methods is presented. 

4.1 Homotopy Methods 

Homotopy methods offer an attractive alternative to more standard approaches of optimal con- 
troller synthesis such as sequential and conjugate gradient methods. The basic philosophy of homotopy 
methods is to deform a problem which is relatively easily solved into the problem for which 
a solution is desired. 

Homotopy (or continuation) methods, arising from algebraic and differential topology, embed 
a given problem in a parameterized family of problems. More specifically, consider sets 0 and Y e 9T' 
and a mapping F : 0 — » Y, where solutions of the problem, 

F(0) = O , (73) 

are^esired with 9 e 0 and F(6) e Y . The homotopy function is defined by the mapping 
//:© x [0, 1] — > such that 


tf(0|.l) = F(0) . 


and there exists a known (or easily calculated) solution, 0 () , such that 


(74) 


h(o q , o) = o 


(75) 


The homotopy function is a continuously differentiable function given by 

H{6{a),a)- 0. ae[0,l] . (76) 

The existence of a continously differentiable homotopy function is assumed, although in many cases 
such as output feedback optimal control design, a continuously differentiable homotopy function 
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may in fact not exist. Also, to be precise, the approach used herein is a continuation algorithm, which 
differs from homotopy algorithms in that the zero curve of a homotopy algorithm is parameterized 
by the arc length. 

Thus the homotopy begins with a simple problem with a known solution, equation (75), which 
is deformed by continuously varying the parameter until the solution of the original problem, equation 
(73), is obtained.-* 51 The power of homotopy methods is that minimization is not strongly dependent on 
starting solution, but depends on local, small variations in the solution. Theoretically, these methods 
are globally convergent for a wide range of complex optimization problems, but in actuality, finite 
wordlength computation often introduces numerical ill-conditioning, resulting in difficulties with con- 
vergence. In light of these numerical limitations, a judicious choice of the initial problem and the associ- 
ated initial stabilizing compensator is necessary for convergence and efficient computation. However, 
the ability to select an initial problem with a simple solution renders homotopy methods more widely 
applicable than sequential or gradient-based methods, which have a stringent requirement for an initial 
stabilizing solution. A new procedure for numerical optimization when the Hessian is ill-conditioned 
or indefinite is presented in appendix A. 

Both discrete and continuous methods are used to solve the homotopy. Discrete methods simply 
partition the interval [0, 1 ] to obtain a finite chain of problems: 

H(0,a ll ) = 0. 0 = a ( ) <a.\ <L <a N = 1 . (77) 

Starting with a known solution at a„.the solution for H(8,a n+ \) is computed by a local iteration 
scheme. Continuous methods involve integration of Davidenko’s differential equation, which is obtained 
by differentiating equation (76) with respect to a, yielding 


d0_ 

da 


' dH\ 1 dH_ 

v d9 J da 


(78) 


Given 0(0) = 6q, this initial value problem may be numerically integrated to obtain the solution atcr = l 
if the solution exists and is uniquely defined. 

In the next section, a continuous homotopy algorithm is presented for fixed-order mixed 7/ 2 / H <*, 
compensator design. The set of homotopy algorithms for synthesizing fixed-order 7/ 2 , H x , and mixed 
H 2 / compensators has been organized into a MATLAB™ toolbox called Fixed-Order Compensation 
of Uncertain Systems (FOCUS). 


4.2 Homotopy Algorithm 

This section describes the continuous homotopy algorithm for mixed // 2 ///„<, control design. 
In essence, a mixed discrete and continuous approach is employed where Davidenko’s differential 



equation, equation (78). is integrated along the homotopy path and at discrete points along the trajectory, 
a local optimization is used to remove integration error. The optimization scheme is a partitioned New- 
ton search method developed for this application and presented in appendix A. Local optimization at 
discrete points along the homotopy trajectory allows a crude integration procedure with large step sizes 
to be employed for efficiently tracking the solution curve. This approach is implemented in the following 
algorithm: 

1. Find initial solution (a=0). 

2. Advance a : 

a \.k =a () + ^ a o,k ■ 


3. Predict 0: 

d (<X\,k) = 0{a 0 ) + Aa 0k O'(a 0 ) . 


where 


da 


c)H \ 1 dH 
dG j da 


4. Check prediction error: 

**(£«) = |/e(0(«l ,k))\ < e 

a. If error less than tolerance, continue. 

b. If not, 0.5A a 0j f — > Aa () . 

c. Increment k and repeat steps 2-A. 

5. Correct with partitioned Newton method to compute local minimum. 

6 . If a = I, stop. Otherwise, go to step 2. 

Various approaches may be taken when selecting the deformation, but the general procedure 
applied in this effort is outlined as follows: 

• Synthesize a low authority H 2 (full order) compensator 

• Reduce the compensator to desired order and transform to canonical form 4 ^ 

• Set y large enough so that the Hi and H M compensators are approximately equivalent 

• Use homotopy to deform the initial low authority, reduced-order // 2 compensator 
into a full authority reduced-order Hi compensator ( Hi homotopy) 

• Deform the full authority Hi compensator into a nearly optimal H compensator 
with 7 approaching its infimum ( H^ homotopy) 

• At discrete values of A, fix y and deform the compensator into the mixed H 2 !H oa 
compensator by varying A (mixed Hi/H^ homotopy). 
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This procedure was chosen because it has been observed numerically that order reduction tech- 
niques tend to work best for low authority LQG controllers. 49 - 52 A canonical form is imposed on the 
compensator structure to minimize the number of free parameters, which in some cases can also lead 
to numerical ill-conditioning. A balancing transformation which does not affect the controller characteris- 
tics relaxes the strict structure in the P° and A 4 ' matrices in equation (29) and improves the conditioning 
of the problem. 


The procedure outlined above separates the compensator synthesis into distinct phases. The initial 
reduced-order full authority compensator is synthesized using the H 2 homotopy, which is then deformed 
into the reduced-order H ^ compensator. During the H x phase, the scalar H 2 norm weight A is fixed 
(as are the plant matrices) and only the norm overbound y is varied. At discrete values, y is fixed 
and A is varied to perform the H 2 norm minimization. Thus, the procedure alternates between the 
and /Gnorm minimization. A similar approach was introduced in reference 15. 

During the homotopy, both the predicted and corrected gains are checked to ensure closed-loop 
stability. After each correction step, the cost gradient is checked to verify descent. During the H x 
homotopy, the solvability of the Riccati equation using predicted or corrected gains must also be checked. 
If any of these conditions are violated, the prediction step size is decreased and the prediction phase is 
repeated. This process continues until the homotopy is completed or until the prediction step size is 
decreased below a prespecified tolerance. 

The following section details the derivations employed in the homotopy algorithm for mixed 
H 2 1 H design. 


4.3 Mixed H 2 IH oo Development 

The homotopy function as well as the gradient and Hessian matrices are determined from the first 
order necessary conditions for an optimal mixed H 2 / //„ compensator given by equations (68)— (72). 
Define 6 be a vector comprised of the free compensator parameters 


6 = vec(G) 


(79) 


where G is the output feedback gain matrix defined in equation (35). The gradient of the cost is 


f(Q , dL ( dL N 

/(0) = — — = vec — — 
dd [dG 


= 0 


(80) 


where dLldG is given by equation (72). 
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The homotopy function is defined as 


H(6.a) 


dL(d.a) 

dd 


= ver 


dL{e,d) 

dG 


= 0 


(81) 


Note that L is now a function of the homotopy parameter a since the system matrices are now functions 
of a . The gradient of the homotopy function is 


V 


H 


(0.a)] = [ 




V H 

v a n 


(82) 


4.3.1 Computation of Hessian 

The derivative of the Nx I vector valued homotopy function. H T (0) = [h l (d),h 2 (6 ), L ,h N {0) 1, 
with respect to the N free parameter vector 6 is the NxN Hessian matrix given by 


— dH 

c)6\ d0 2 dQ]y 


where 


c)H 

ddj 



and using equation (72). 


dH 

ddj 


vec 


yddj 


: {2[(d,^D 12 GC 2 -D&Ci-Bj qJ)LcJ 


+ (tof p D ip GC 2 - ADfpCp - bJ L p)xcJ ] } 


= vec 


A2A 2 g U) c 2 -bJq { J ) 



+ {px 2 D n GC2 - A 2 A - b[ Q^L^cT 

+ {XDf p D X pG {j) C 2 - bJ l\ i^xcj 


+ (to l T p D [p GC 2 -tof p C p 



(83) 


(84) 


(85) 


(86) 
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The derivatives with respect to 0 are denoted 


/*)(./) = itl 

y ’ d0 ; 


For 0j = t > ik , 


G \.n = z^ = E 

d6j 


c)G 


which is a matrix of zeros except for a one in the { i , k) element. 

To obtain expressioi 
with respect to 6 . to obtain 


To obtain expressions for l} J \ Qli\ I^j}\ and X^\ differentiate equations (68 )— ( 7 1 ) 


0 = (i4 + y- 2 BB T QjjL U) + 6 j) [A + y~ 2 BB T Q^ 


(87) 


( 88 ) 


+ 


A {,) + y~ 2 BB T Q { J ] V + l[ A + y~ 2 BB T Q ( J^ T 


(89) 


o = (ii + Y~ 2 BB T Qjf Q { J } + Q { J\a + Y~ 2 BB T Q^) 


+ 


A U)T Q^+Q^A (i K[c r cf ] 


(90) 


0 = A T L ( f + 6pA + 


A U)T L p + L p A U) + X(clc p f j) 


(91) 


0 = AX (/) + X^A T + 


A^X + XA^ 7 + [b p bJ } j j) 


(92) 


Derivatives of the closed-loop matrices are obtained from equations (64)-(66) and are given by 




(93) 


(bb t ) {j) = o, (b p b 7 ) u, = 0 




(94) 
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T 


C T C 



[c 2 G U)T D l T 2 D l 2 GC 2 ^j + [c 2 G ( ' j)T D^ 2 D l 2 GC 2 ) 
-[c^D n G {i) C 2 ] - [cf D V 2 G {J) C 2 j 7 


( 95 ) 




-0 lp G ( y) C 2 


C /? -D 1/ ,GC 2 )+(c / ,-Z) I/ ,GC : x7 ' 


-D| /; G (/) C 2 


(96) 


4.3.2 Computation of H a 

Similarly, the derivative of the homotopy function with respect to the homotopy parameter, a, 
is the iVxl vector 


V„H r = dH 


da 


(97) 


and 



(98) 


d „ 

reel —^7 

da 


D ' 2 D\ 2 GC 2 - D\ 2 C] - B 2 QjjLC 


t?T 

1 


+(U>f p D lp GC 2 


folpCp 


bJ l p )xc t 


]}) 


(99) 


= vecl 2 


j D 12 D 12 GC 2 + d[ 2 D 12 GC 2 - D t2 C| - bJqJ\lc{ 


+(d 12 d 12 gc 2 - d 12 q - bJq^lcT 
+(Xd{ p d Xp gc 2 - XD? p c p - bT L p 

+A(B,;D Ip CC 2 + D^ApGC, -5,^c p ))xc% r ]| , 
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( 100 ) 



where 



da 


( 101 ) 


Implicit in these equations is the assumption that only the system matrices D (2 and D \p an d the 
parameters A and y are functions of a. In general, the homotopy can be performed with any/all system 
matrices deformed. 

The derivative terms in equation ( 100) depend on the deformation undertaken in the specified 
problem, i.e., the initial and final problem. Suppose that the deformation of the matrix Z)p is prescibed 
to be 


D\l( a ) ~ ^I2.()( a ) + a (^\2,f( a )~ ^12.()( a )) ' (102) 

where the 0 and/ subscripts indicate the initial and final system matrices, respectively. It follows that 

D^i — D\2'f ~ 7 / 2.0 • ( 103 ) 

The derivatives of other plant matrices and the parameters A and y are determined accordingly. 

To obtain expressions for the derivatives of L, Q x , L p , and X with respect to a, equations (68)— (7 1 ) 


are differentiated, resulting in the following: 

T 

0 = [A + y- 2 BB r Q^)L + L[A + Y~ 2 BB T Q oo ) ) +(rL + LT r ) (104) 

0 = (4 + y~ 2 BB T Qjf a 00 + Qoo(A + y~ 2 BB T Qj) -{ly^yQ^BB 7 Q x ) ( 1 05) 

0 = A 7 L p + LpA + ^ A 7 L p + L p A + AC 7 Cp + AC 7 C p + AC 7 C p j (106) 

0 = AX + XA T +^AX + XA T +B p Bj y +BpBl^j , (107) 

where 

T = A- 2 y~ 3 y BB 7 Q x + Y~ 2 (bB T + BB 7 Q x + BB 7 Q x j , (108) 
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and from equations (64)-(66) 


4 - 4 - B~,GC-> - /m;c\ 


(109) 


B P = Bp 


5 = 5, 


Cp -C p -D lp GC 2 -D Xp GC 2 


C = C] - D X2 GC 2 -D\ 2 GC 2 , 


where from equations (58 )— (6 1 ), the augmented matrix derivatives reduce to 


A/ - /V) 0 

-/V (C 2 ,/' o) 0 


B p,f ~ B p . o 


P ~ N °{ D 2p.f -D 2 p,()) 


B l.f~ B l,0 

-a/ () (d 21 / -d 2 io) 


Bif — Z?2,() 




^/» = [ C /». / _C 5.0 °] 


(118) 


Q -[Cl,/ _ Q,o °] 


(119) 


C 2 =[0 0] 


( 120 ) 


A p - [a p, 


1 P.f ~ "lp.0 


( 121 ) 
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( 122 ) 


D \2 ~[ D \2,f - °12.0 °] • 


The presence of the zero subblocks significantly enhances the computational efficiency of this 
approach. When implementing the procedure described at the beginning of this section, the above equa- 
tions may be further specialized. In the initial H 2 homotopy procedure, the initial and final plant matri- 
ces are the same and the homotopy is performed only on the measurement and process noise intensities, 
D ]2 and £> 21 - Hence A, Bi, C|. and Ci are identically zero. For the mixed Hj ! homotopy, the H M 
and Hj homotopies are performed distinctly, which simplifies the computations significantly since the 
plant matrices remain fixed and only 7 or A are varied at one time. 
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5. DESIGN EXAMPLES 


5.1 Example 1: Flexible Four-Disk Example 

To demonstrate the homotopy algorithm applied to optimal controller synthesis, the four-disk 
example originally described in reference 53 and more recently by numerous others 49 will be used. 

The four-disk model used in the example problem was derived from a laboratory experiment and repre- 
sents an apparatus developed for testing of pointing control systems for FSS’s with noncolocated sensors 
and actuators. Four disks are rigidly attached to a flexible axial shaft with control torque applied to 
selected disks and the angular displacement of selected disks measured. The plant parameters are taken 
from reference 1 . 

The seminal paper dealing with mixed H 2 / H ' x design addressed the case where vc = w p in 
equation (58) with results from the four-disk problem given for the full order case.l In this section, the 
one input, two output case will be repeated with the FOCUS algorithm as well as a two input, two output 
case with full- and fixed-order compensators to demonstrate the capabilities of the homotopy algorithm. 

Table 1 presents a comparison of the results from the FOCUS algorithm and the results published 
in reference 1. In table 1, the "BH" subscript indicates results from reference 1 and the “F” subscript 
indicates results from FOCUS. Gaps in the columns of table 1 denoted “BH” correspond to values of y 
for which results were not published. The absence of data generated by FOCUS for y = 2 and 1 .5 is due 
to the fact that y is larger than the maximum oo norm as described in the next section. A key distinction 
between the formulations of the mixed H 2 / optimization problem used herein and that of 

reference I is that while their formulation minimizes an overbound on the // 2 norm, the formulation 
utilized in FOCUS minimizes the actual // 2 norm. Consequently, as shown in figure 4, the FOCUS 
algorithm results in smaller // 2 norms for a given y with the gap in the H x norm overbound smaller 
than the reference 1 results in the meaningful region for y, as described in the next section. Whereas 
7 = 0.52 was the smallest 7 value reported in reference I, 7 values of 0.5, 0.4, and 0.3 are utilized with 
FOCUS, as indicated in table 1 . A quasi-Newton technique based on the Broyden-Fletcher-Goldfarb- 
Shanno (BFGS) variable metric algorithm was used to extend these results for a full-order mixed 
// 2 / Woo controller, achieving 7 = 0.29, corresponding to an actual H ^ norm of 0.28. 54 

This example also demonstrates the synthesis of fixed-order mixed controllers with singular 
constraints. In practice, compensators representing the extreme values of 7 or A will not typically be 
used since either performance or robustness would be severely diminished. A compromise value in the 
“elbow” of figure 4 would typically be chosen. While the FOCUS results in table 1 are not significantly 
better than the results of reference I , the example demonstrates that the FOCUS code performs satisfac- 
torily and lessens the gap between the overbound and the H ^ norm. 
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Table 1 . Comparison of FOCUS and reference 1 results. 
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T Z pwp 2 BH 

T 

zpwp 2 F 

2.0 

1.18 


0.382 


1.5 

1.06 


0.389 


1.0 

0.855 

0.947 

0.410 

0.398 

0.9 

0.797 

0.889 

0.420 

0.403 

0.8 

0.732 

0.776 

0.432 

0.421 

0.7 

0.661 

0.668 

0.451 

0.438 

0.6 


0.589 


0.460 

0.52 

0.511 

0.520 

0.512 

0.508 

0.5 


0.5 


0.518 

0.4 


0.398 


0.572 

0.3 


0.300 


0.799 



Figure 4. H x versus cost for four-disk example. 


This example serves to illustrate some interesting features of the mixed Hj / formulation 
implemented in FOCUS as well as distinctions from other formulations. The formulation implemented 
in FOCUS is a method for generating suboptimal H x controllers of fixed order that are subject to an H 2 
constraint. The cost functional for the mixed H 2 / H problem can be written as 

J mix = J 00 +XJ 2 * 023) 
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where 


./oo = tracejeA^fifi 7 } 

(124) 

J 2 = trace jxC^C^J , 

(125) 


subject to the corresponding Riccati and Lyapunov equations. The resulting lagrangian is given by 
equation (67). 

Minimization of J a 0 results in an H M controller with an upper bound on the oo norm given by 
the oo norm of the Hi controller (as y — » », the H 2 compensator for T. w is recovered by minimizing 
./,*,). By successively lowering gamma, the minimum norm controller for T :w is obtained. Minimi- 
zation of J 2 results in the optimal Hi compensator for T :pwp . So when nonzero A is used in J nnx , the 
H x cost functional imposes an additional constraint on the H^ norm of T . w , and for large y. the 
necessary conditions for the mixed problem yield the simultaneous solution of two Hi problems. By 
increasing the Hi weight. A , for a fixed y, is reduced while ||T- m ,|| oo approaches the gamma 

overbound. At that point, the minimum Hi controller for T :pwp such that ||7’- v ,.|| oo < y is obtained. 

For this example problem, the optimal Hi controller for T. pwp results in |7\. ,,1 = 1.392 and 
||r rvi .||^ = 0.3786. y values larger than 1 .392 are not meaningful for our formulation because the optimal 
Hi compensator for T zpwp satisfies the y constraint on ||7'- H ,|| (x) . In that case, the H^ norm constraint is 
inactive in the mixed-norm optimization. However, since the reference 1 formulation seeks to minimize 
the Hi norm overbound, it is possible for their formulation to generate meaningful solutions for 
y > 1 .392. When using FOCUS, the first step should be to establish an upper bound on y by computing 
|77 m .L resulting from the optimal Hi compensator for T. pwp . 

5.2 Example 2: Building Control Example 

An interesting example problem is that of controlling the vibration of a building subject to an 
earthquake excitation. The problem from a controls point of view is the need to develop a controller that 
can reliably accommodate the uncertainty in excitation that is characteristic of earthquakes, while at the 
same time handle the presence of uncertainties caused by inelastic structural response. This example 
examines design approaches which achieve nominal performance only ( H 2 ), robust performance 
{/J -synthesis), and nominal performance/robust stability (mixed H 2 / jtf), applied to the problem of 
building structural control. The challenge is to achieve the highest attainable level of RMS performance 
for a specified bounded set of uncertainties. A comparison of these controller design techniques is given 
based on an experimental model of a three-story tendon controlled structure at the National Center for 
Earthquake Engineering Research. 55 The laboratory structure is depicted in figure 5. 
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Figure 5. Structural control experiment. 


5.2.1 Control Design Model 

The evaluation model for performance assessment is a 20-state model from reference 55. ob- 
tained by system identification experiments on the laboratory structure. A 6-state nominal design model 
was obtained by balancing and residualizing the 20-state evaluation model, retaining modes at 2.268. 
7.332, and 12.240 Hz. Inputs to the model consist of the ground acceleration disturbance, ,v y , sensor 
noise, and the tendon control input, u. Performance outputs include the weighted displacement of the 
three floors relative to the ground, zp , and the weighted control force, zu. The measurement output, y, is 
the absolute acceleration of each of the three floors. All units are in volts. Figures 6 and 7 illustrate the 
fidelity of the reduced-order design model by comparing the frequency response of the design model and 
evaluation model and are typical of the other input-output pairs. Figure 6 presents the transfer function 
from the control input to the relative displacement of the third floor. The frequencies of the first three 
dominant modes are matched well, although some error exists with the modal gains. Figure 7 shows 
excellent matching with the transfer function from the disturbance input to the relative displacement of 
the third floor. 

To provide robustness to model errors, the design model is extended to include parametric uncer- 
tainty within the control bandwidth in the form of errors in the modal damping and frequency squared 
terms as introduced in reference 17. Uncertainty will only be used for the natural frequency squared 
terms in this design, but for completeness, the formulation for uncertain modal damping will also be 
presented. Although the uncertain natural frequency square terms are real parameters, using a complex 
uncertainty also accounts for variations in modal damping if a hysteretic damping model is assumed. 

In modal form, the nominal A matrix for a second order system is written 


A o 


0 1 

-o) 2 -2 C,a> 


(126) 
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Figure 6. Comparison of design and evaluation models for relative 
displacement of floor 3 to control input. 



Figure 7. Comparison of design and evaluation models for relative 
displacement of floor 1 to disturbance input. 


Introducing multiplicative uncertainty in the modal frequency square and modal damping terms 
results in 


A = 


0 

~(o 2 (\ + <5i) 


1 

-2£(o{\ + S 2 ) 


(127) 


= Aq + AA 


(128) 
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where 


AA = 



0 

-2i^cod 2 


(129) 



' 0 “ 

-co 2 0 

+ <5~> 

' 0 ' 

<5, 




1 



_ 1 


-2 C«] • 


For a system with n total inodes and ni uncertain inodes, A = /)<) + X/=i A A,-, and 


( 130 ) 


M = { e 2i)8li(- (0 t){ e 2i-lj r + ( e 2i)^2i{~^i (0 i)( e 2i) 1 


(131) 


where (e ; j is the jth standard basis vector for g^ 2 ". Defining A to be the set of indices of uncertain modes 
allows the plant with uncertain natural frequency square and damping terms to be as shown 
in figure 8 with the following definitions: 


A A lw - A A ld - E 2k , A A rw - AA rd ~ -DE1E 2 \ 

Q = diagf(0^,j , D = diag[2^ ( jJ, V/'=I,2,K ,m , 


= e 2 


2k -K2A(1) c 2k(2) 


^ ^ A /- 


k(m 


.]• 


(132) 

(133) 

(134) 



Figure 8. Plant with uncertain modal damping and frequency square terms. 


In addition to the uncertainty in the modal frequency square terms, an additive uncertainty 
is included to represent model error outside the control bandwidth. This type uncertainty model forces 
the controller to gain stabilize the high-frequency modes that were truncated from the evaluation model. 
For robust control design, the baseline uncertainty model included 5-percent uncertainty for the natural 
frequency square error {Wm = VO. 05 ) and an additive uncertainty weighting function given by 


W, 


additive 


= 6.4 


(* + 5) 3 
(s + 200) 3 


K add * ^add 


( 135 ) 
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Figure 9 presents the transfer function from control input to the acceleration of floor I and the corre- 
sponding additive uncertainty weight. In order to balance the plant for improved numerical results, 
the additive uncertainty model is realized as the frequency-dependent term. W a(](l , and the constant 
gain term. K (l(h/ . 



Additional inputs for the robust control design generalized plant include inputs associated with 
the additive uncertainty, vr r and the modal frequency uncertainty, m’ wj . Additional outputs include those 
associated with the additive uncertainty. z (/ . and the modal frequency uncertainty, z . The uncertainty 
block has the structure 


< 5 , 


A = 



(136) 


with A 4 g C‘ 4aI and A^ g C 4a4 . Figure 10 illustrates the generalized plant for robust control design. 

5.2.2 Control Design Results 

This section presents results of the design approaches for nominal performance ( // 2 ), robust 
performance (/t -synthesis), and nominal performance/robust stability (mixed // 2 //t) for the building 
structural control problem. For evaluating the nominal performance of these designs, performance is 
defined by Vz, the trace of the output covariance where the outputs are the three relative floor displace- 
ments, and Vn , the trace of the control input covariance. 
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Figure 10. Generalized plant for robust control design: Building control example. 


For the //? nominal performance design, the disturbance input and performance output 
vectors are 


noise 


- 

~ _ 




- 

L J 


~n 


(137) 


with the control weight. Wu = ^Jp . the weight on the relative displacement of each floor, Wp = 25, 
the sensor noise intensity, Kn - 0.001 , and the intensity of the ground disturbance, Kd = 0.001 7. Kd 
was chosen to match the dc intensity of a reference earthquake excitation known as the Kanai-Tajimi 
(K-T) spectrum. 55 Control authority was varied in the design process using the scalar p . 

For the p -synthesis design, the corresponding disturbance and performance vectors are 


w n, 


“ m 

w a 


-a 

noise 


~ a 

h 


i 

1 1 

1 


Aset of p controllers of varying control authority was designed by fixing Wp and varying p . At each 
level of control authority, the relative weighting between z p and z u was determined by p and with the 
uncertainty weights fixed, z= z p , z„ was scaled (corresponding to the performance block for the p 
design) to achieve a p measure of one. By scaling the performance variables in this manner, at each 
level of control authority the controller was designed to maximize performance while providing a fixed 
level of robustness. Hence, a consistent comparison of control approaches was made from a robustness 
perspective. First order D-scales were used for each p controller design, resulting in p controllers with 
19 states computed using the MATLAB™ p -Analysis and Synthesis toolbox. 56 
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Finally, a set of mixed H 2 I p controllers were designed with fixed controller dimension of sixth 
order using the homotopy algorithm presented in section 4. In order to trade between nominal perfor- 
mance and robust stability, the H 2 subproblem is defined for nominal performance as above and the 
p subproblem accounts for the additive and parametric uncertainty models. This paradigm is an effective 
means of exploiting the inherent trade between robustness and performance in mixed-norm design to 
separate the two competing design objectives. The H 2 subproblem is defined by the performance vari- 
ables and the p subproblem is defined by the variables associated with the uncertainty structure. The 
subproblems are defined by the inputs and outputs 


’ w w " 

r l - 


VI ’7 — 

noise 

z 2 — 

Z P 

. w a _ 

~a 


. *8 . 


_ 


and the D-scales for the p subproblem are obtained from D-K iterations for T- |M| . 

Figure 1 1 presents the mean-square (MS) nominal performance curves for each control design 
method. The robust control designs are for the baseline uncertainty model (which has 5-percent uncer- 
tainty in the natural frequency square parameters and the additive uncertainty). The costs are computed 
for the closed-loop with input noise filtered through the K-T spectrum. H 2 design costs are computed 
for both the design and evaluation models to illustrate the limitation on achievable performance due to 
model error. Although the cost curve evaluated with the design model extends to high control authority 
levels, the maximum performance with the evaluation model is obtained at p - 15.63 (indicated by “o” 
in fig. 1 1 ). The loop closed with the H 2 controllers and the evaluation model are unstable for smaller 
values of p. This cost comparison also indicates that for control authority levels lower than the instabil- 
ity level, the actual performance is almost identical to the design model performance. 

Figure 1 1 also indicates the loss of MS performance that is incurred in exchange for robust 
performance. As a basis for comparison, the set of p designs is evaluated in terms of MS performance. 
A substantial gap in performance exists between the H 2 and p designs since the p designs achieve 
a given level of output performance at a higher control cost than the H 2 designs. However, the mixed 
H 2 / p designs effectively recover the MS performance of the H 2 designs while providing the same 
level of robust stability as the p designs. The mixed H 2 Ip design procedure provides performance 
comparable to H 2 design while overcoming the major shortcoming of H 2 design, namely a lack of 
stability robustness. 

The impact of uncertainty on performance in the mixed-norm design setting is evident in 
figure 12 where a set of mixed-norm designs for 10- and 20-percent parametric uncertainty are evaluated 
in addition to the baseline 5-percent parametric uncertainty mixed-norm designs. As the level of robust- 
ness increases, performance is sacrificed as indicated by the upward shift in the performance curve. 

A cursory comparison of figures 1 1 and 12 indicates that the mixed-norm controllers designed for 
10- and 20-percent parametric uncertainty yield comparable performance to the p controllers designed 
for 5-percent uncertainty. Hence, the mixed-norm designs provide more robust stability for a given level 
of performance than the p controllers. Note that these comparisons are for nominal performance and 


38 



may not hold for robust performance. In these analyses, the additive uncertainty is held fixed since it is 
defined with respect to the model and serves only to force the controller to roll off and gain-stabilize the 
high-frequency unmodeled dynamics. 



Figure 1 1 . Mean-square performance comparisons. 



Figure 12. Impact of uncertainty level on mean-square performance. 
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Robust stability of each design is evaluated using mixed /A analysis where the parametric uncer- 
tainty is considered real and the additive uncertainty complex. As a result, the mixed /A measure is a less 
conservative measure of robust stability. Figure 13 plots the ja measure for the set of H 2 controllers for 
varying authority levels as a function of parametric uncertainty level. This plot should be interpreted 
as indicating the magnitude of perturbation required to destabilize the closed-loop. From equation (19), 
a jA measure < 1 indicates robust stability is guaranteed for all plants in the uncertain set. For a control- 
ler associated with a [A measure of / 3 . the system will be unstable for lA^ >j^. The H~, designs are 
robust with respect to the uncertainty model only for very low authority controllers. Figure 13 illustrates 
the well known property of H 2 controllers that as control authority increases, the sensitivity (in terms 
of stability) to model error increases. This figure also indicates that the jA measure is relatively insensi- 
tive to different levels of parametric uncertainty at high control authority levels which indicates that the 
additive uncertainty dominates the stability analysis. Only at low authority levels are the H 2 designs 
sensitive to parametric uncertainty. Since control bandwidth is proportional to the authority level for 
these // 2 designs, the higher authority controllers interact with and destabilize the unmodeled modes. 



Figure 13. Robust stability analysis of H 2 controllers. 

Robust stability analyses of the mixed-norm designs for 5-percent, 10-percent, and 20-percent 
parametric uncertainty are shown in figures 14-16. For the 5-percent uncertainty design, an H x over- 
bound of one was achieved. Although robust stability is not guaranteed for levels of uncertainty above 
5 percent, the /A measure for 20-percent parametric uncertainty is <2, which is roughly 3 times better 
than the H 2 designs. It is also interesting to note that the /A measure for the mixed-norm design 
is sensitive to differences in parametric uncertainty and is relatively insensitive to the control authority, 
which is opposite the characteristic of the H 2 designs. As a matter of fact, the (A measure decreases 
slightly with control authority for the mixed-norm designs. Somewhat different behavior is observed 
with the mixed-norm designs for 10- and 20-percent parametric uncertainty. The mixed-norm design 
set for 10-percent parametric uncertainty used an overbound of 1.3, so robust stability is not fully 
guaranteed for 10-percent variations in the uncertain natural frequency. The peak /a measure in figure 15 
is 1.26. Similarly for the mixed-norm design with 20-percent parametric uncertainty, an H x overbound 
of 2. 1 was used and the peak fA measure is 1 .75. These two designs have a characteristic behavior more 
similar to the H 2 designs in that the fA measure is more sensitive to control authority than parametric 
uncertainty level. However, the variation with control authority is significantly less than the H 2 designs. 
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Figure 14. Robust stability analysis of mixed-norm designs performed for 5-percent uncertainty. 



Figure 15. Robust stability analysis of mixed-norm designs performed for 10-percent uncertainty. 



Figure 16. Robust stability analysis of mixed-norm designs performed for 20-percent uncertainty. 
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Fora second-order system with an uncertain natural frequency square parameter (equation (127) 
with S 2 = 0), if <5] is considered a real parameter, the uncertain system will be stable when c>i > -I . 
However, if <5[ is a complex variation, the system is stable only when |5]| < 2£.- £ ' 7 Thus representing the 
real parameter uncertainty as a complex variation introduces significant conservatism in the control 
design. The impact of this is evident in the mixed-norm design with 10- and 20-percent parametric 
uncertainty. The homotopy began with a fixed-order p design for 71 M . which exists because of the 
artificial destabilizing effect of the complex parametric uncertainty. For the 10-percent uncertainty level, 
the fixed-order H design with the D-scales from the full-order p design resulted in a minimum H x 
norm of 1.2539. For 20-percent uncertainty, the minimum norm is 2.0463. Since these designs are 
for the H rx subproblem only, they represent a lower limit on the // M norm for the mixed H-> / p designs 
and are used as the initial points for the A homotopies in the mixed-norm designs. 

The complex p measure of each mixed-norm boundary controller is only very slightly less than 
the H x norm overbound, indicating that the D-scales obtained from the full order D-K iteration for 
T. Wl for each uncertainty level are nearly optimal for the sixth order mixed-norm controllers. Had this 
not been the case, the D-scales could have been optimized for a mixed-norm boundary controller, fol- 
lowed by another fixed order controller optimization step. 

This design example also serves to illustrate the characteristics of the homotopy algorithm along 
the solution path. The mixed H 2 / p design begins with a full-order p design for 71 M . followed by 
order reduction and a y homotopy to obtain the minimum fixed-order p design for 71 M . . This corre- 
sponds to a mixed H 2 / fA design for A = 0, from which the A homotopy begins by fixing the y 
overbound slightly higher than the minimum y and incrementally increasing A. A is increased until 
the H 2 norm of 71 H ., reaches a minimum and the H x norm of 71 H . is approximately equal to the 
yoverbound. At that point, one boundary design is obtained for a given value of p . The set of boundary 
designs are computed by varying p from this point. 

As the A homotopy progresses, the H 2 norm decreases from the minimum fixed-order 
fA controller value to the minimum H 2 norm subject to the H ^ norm overbound. The variation of the 
H 2 norm during the A homotopy for the 5-percent uncertainty level is shown in figure 17 along with 
the Il 2 norms obtained with the fixed-order p design and the full-order H 2 design. The mixed-norm 
design reduces the H 2 norm of the fixed order p design by 50 percent while providing the same level 
of robust stability. Conversely, the mixed-norm design only increases the H 2 norm by 10 percent over 
the full-order H 2 design in exchange for robust stability. Note also that the most significant H 2 cost 
reduction occurs in the first few iterations which indicates that the entire A homotopy does not have to 
be performed to make substantial improvements in H 2 cost. Figure 18 shows the constraint on the H x 
norm is satisfied as the H 2 norm decreases to the minimum. 
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The spikes in figure 1 8 are artifacts of the numerical optimization algorithm and are indicative 
of the numerical sensitivity of the homotopy algorithm. The controller canonical form tends to be poorer 
conditioned than a nonminimal realization which makes optimization of the fixed-order controller 
parameters difficult. After order reduction and transformation of the controller to canonical form, the 
initial gains are typically not optimum and the Hessian is often ill-conditioned and indefinite. Standard 
Newton optimization methods fail when the Hessian is ill-conditioned and indefinite. When using 
a standard Newton optimization method, the homotopy would terminate prematurely due to the ill- 
conditioned and indefinite Hessian. This numerical sensitivity motivated the development of a numerical 
optimization method that accounts tor ill-conditioned and indefinite Hessians. Implementation of this 
optimization algorithm, described in appendix A, resulted in a numerically robust algorithm which 
converged along the homotopy path. 

As for the homotopy along the boundary for variation in p, a typical gain trajectory is shown 
in figure 19 which is the second element of the parameter vector corresponding to the ( 1 ,2) element 
of the gain matrix G. Note that the gain variation is not monotonic. The trajectory of control gains 
in essence provides the mechanism by which the gains of a multivariable controller may be “dialed in.” 
much the same manner as the dc gain of a classical controller. This enables an incremental implementa- 
tion which may be monitored to determine the onset of instability and the maximum achievable nominal 
performance for a given design model. Hence, this example illustrates how the mixed-norm design 
procedure may effectively be used as a means of tuning multivariable controllers in orbit to achieve 
maximum performance. 



P 

Figure 19. Gain variation along boundary homotopy with 5-percent uncertainty. 
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5.3 Example 3: Flexible Space Structure Experimental Example 


5.3.1 Facility Description 

The Controls/Structure Interaction Ground Test Facility (CSI GTF) at the NASA Marshall Space 
Flight Center (MSFC) has been developed for experimental investigation into the control, system identi- 
fication, and dynamics of an FSS. Characteristics of FSS’s which make modeling, simulation, and 
control design challenging are embodied in the facility. The experiment is a large, flexible structure 
which has numerous low-frequency, coupled, tightly spaced, lightly damped modes. The present con- 
figuration includes a flexible, deployable boom which once flew on the Shuttle for the Solar Array Flight 
Experiment (SAFE). Initially, the facility was designed as a GTF for the Control, Astrophysics, and 
Structures Experiment in Space (CASES) flight experiment. 

As shown in figure 20, the test article is inverted and vertically suspended from a platform 132 ft 
above ground level. A disturbance system provides two translational degrees of freedom to the base of 
the structure (the base refers to the portion of the experiment at the top platform). A simulated mission 
peculiar equipment support structure (MPESS) interfaces the disturbance system with the test article 
to simulate a flight experiment interface between the Shuttle, MPESS, and the payload. The test article 
consists of a 105-ft boom which supports a simulated occulting plate at the boom tip. The control objec- 
tive of the flight experiment is to maintain alignment of the tip plate with a detector at the base of the 
boom on the Shuttle; this would allow the occulting plate to point towards a star to perform an x-ray 
experiment. Similarly, the ground experiment strives to maintain alignment of the tip plate with the 
simulated detector at the MPESS. 

Control authority is provided by angular momentum exchange devices (AMED's) and bi- 
directional linear thrusters (BLT's). Each AMED package consists of two motors attached to reaction 
wheels, two two-axis gyros, and associated electronics. One AMED package is located at midboom 
and a second AMED which is augmented with a third reaction wheel is located at the tip. Each AMED 
motor has a peak rated torque of 290 oz-in. Two single-axis bidirectional linear thrusters are located 
at the boom tip which have a force capability of ± 2 lb up to 10 Hz. The shakers are considered actuators 
for disturbance generation and not for control. The measurement system includes sensors available for 
feedback control, disturbance, and safety monitoring. Control sensors include three single-axis acceler- 
ometers at the base (MPESS), two dual-axis gyroscopes (one redundant axis) in the midlength and tip 
AMED packages, three single-axis tip accelerometers, and a three-degree-of-freedom tip displacement 
sensor (TDS). Hence there are seven control actuators, 15 control sensors, and two programmable 
disturbance inputs. The CSI GTF computer system consists of a Sun host and a separate real-time 
system. Real-time control is provided with the capability of 64 sensor measurements, 64 actuator com- 
mands, and up to a 100th order controller at a rate of 250 Hz. 

5.3.2 System Identification 

The CSI GTF exhibits the pathologies that make control of FSS’s challenging. System identifica- 
tion and control design for this structure is difficult due to the modal density and light damping at low 
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Figure 20. CSI ground test facility. 


frequencies. Since the tip displacement is dominated by the low frequency bending modes, the model 
fidelity must be high in this regime. Modal density at low frequency also leads to high-order design 
plants. 

The Eigensystem Realization Algorithm with Data Correlation (ERA/DC) 58 was used to obtain 
control design models for the CSI GTF. Time response data are used by ERA/DC to generate a discrete- 
time state space realization of the system. Inputs for system identification were the two disturbance 
shakers, the two tip thrusters, and the tip z-axis reaction wheel. In the initial iteration of system identifi- 
cation, each actuator was individually excited with an 80-sec burst of uncorrelated random noise with 
impulses interspersed every 20 sec. Comparison of the resulting models with frequency responses of the 
actual data indicated that several significant modes below 1 Hz were not accurately identified. Impulsive 
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inputs primarily excited the dominant first bending modes. To provide more energy at the frequencies 
where the unidentified modes were expected (from modal testing), a set of inputs were formed by sum- 
ming sine waves at the 10 discrete target frequencies with a random phase. An uncorrelated noise se- 
quence filtered with two poles at 25 Hz was added to each input in an effort to better facilitate identifica- 
tion of zeros. For one experiment, each actuator was excited simultaneously for 200 sec and the response 
was measured at the feedback control sensors. This procedure was repeated four times. The 200-sec 
response time histories were sampled at 250 Hz and postprocessed by removing the mean, low-pass 
filtering with four poles at 5 Hz and decimating to a 10-Hz sample rate. 

To generate control design models from the time history data using ERA/DC. the System/ 
Observer/Controller Identification Toolbox for MATLAB™ was used.-"' 9 A 100-state model was obtained 
which identified the significant modes in the frequency range of interest for control design. This 
lOOth-order model was reduced to a 40-state model by transforming to a balanced realization and trun- 
cating the states with low controllability/observability grammians. The modes above 3 Hz were then 
residualized resulting in a 26-state nominal model of the plant dynamics used for control design. Table 2 
presents a comparison of the modes below 3 Hz from the 40-state system identification model with those 
generated by modal testing 6() and finite element modeling^ and a description of each mode. Note that 
the modal test and system identification tests used different excitation and sensing. 

As table 2 indicates, the structure exhibits a high degree of modal density with closely spaced, 
lightly damped modes, which taxes the system identification procedure. However, good results were 
obtained by ERA/DC and the resulting design model accurately describes the dynamics of the CSI GTF 
in the desired frequency range of <3 Hz. Due to the offset of the center of mass from the boom axis in 
the .v direction (as well as coupling with the tip suspension system), bending in y couples strongly with 
torsional motion. Accurately representing this dynamic coupling in the identified model is critical for 
control design and is an indication of the fidelity of the ERA/DC identified model. The third and fourth 
.v and v bending modes involve coupling with tip plate bending also. 


Table 2. Comparison of experimental and analytical frequencies. 


System ID 

Modal Test 
Frequency 
(Hz) 

Finite Element 
Frequency 
(Hz) 

Mode 

Description 

Frequency 

(Hz) 

Damping 

<%) 

0.118 

0.506 

0.112 

0.08 

y 1st Bending 

0.119 

6.48 

0.120 

0.09 

x 1st Bending 

0.215 

0,904 

0.210 

0.16 

1st Torsion 

0.325 

57.3 




0.535 

1.07 

0.520 

0.56 

x 2nd Bending 

0.554 

0.604 

0.530 

0.58 

y 2nd Bending 

0.953 

3.60 




1.214 

38.8 

1.391 

1.23 

x 3rd Bending 

1.539 

17.7 




1771 

3.15 




1.872 

2.76 

1.868 

1.82 

y 3rd Bending 

2.061 

13.2 




2.839 

9.45 

2.802 

2.73 

y 4th Bending 

3.073 

3.58 

2.995 

3.28 

x 4th Bending 
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5.3.3 Control Design 

The control objective for CSI GTF is to maintain alignment of the tip plate with respect to a 
sensor at the base of the boom when subjected to a disturbance force (DSy) applied at the base. To 
accomplish this objective, the measurements used for control are displacement sensed by the tip dis- 
placement sensors ( TDS.x and TDSy) and the three rotational rates sensed by tip rate gyroscopes ( TGx , 
TGy, and TGz). Sensor noise is added to the plant output for control design. The BLT and the tip r axis 
reaction wheel (MC5) comprise the set of actuators used for feedback control. Due to the large general- 
ized mass of the first x and y bending modes, the BLT’s are the only actuators that have sufficient control 
authority to affect these modes. 

The design results which follow are based on the mixed H 2 / control methodology to achieve 
high nominal performance as measured by the MS alignment error of the tip plate while providing 
stability robustness in the presence of model errors. A compensator dimension of five was chosen for 
these design comparisons. Enforcing a constraint on the compensator order of five states is an extreme 
case chosen to highlight the effect of compensator dimension for fixed- and reduced-order designs. Five 
states result in a true canonical form for the compensator (the plant has five measured outputs) and also 
reduces computation. A major factor contributing to the computational burden of the mixed-norm 
homotopy algorithm is computing the Hessian. Since the Hessian is computed one row at a time and the 
number of rows is a multiple of the compensator dimension, keeping the compensator dimension small 
is important for computational efficiency. 

A set of full-order H 2 designs provides a baseline for nominal performance evaluation with 
weighted control and weighted tip displacements as performance outputs and disturbance force and 
sensor noise as inputs. With reference to figure 21. the H 2 subproblem is defined by 


U’2 — 

noise 

~2 — 

Z P 


DSy 


_ 


(140) 


The baseline nominal performance design uses the weights KDSy = 25, Wp = 100, Kn = 0.001 , and 

= yp where p is scaled to vary control authority. Figure 22 shows the performance of the full-order 
baseline H 2 designs and reduced fifth order H 2 compensators evaluated by comparing output and con- 
trol cost with p varied from 1 ,000 to 100 (p is inversely proportional to control authority). The output 
and input costs plotted in figure 22 are the traces of the output and control covariances. Note that the 
solid line corresponding to the fifth order H 2 controllers obtained by balanced model reduction indicates 
that the order reduction process does not preserve closed-loop stability with the design model for higher 
levels of control authority (corresponding to p values <600). 
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Phase I: 


In the initial mixed-norm designs, the H ^ subproblem employed a frequency-dependent additive 
uncertainty, W (uj(jir nr , to constrain control bandwidth and a constant multiplicative uncertainty at the 
plant input, W mu / t , to account for model error in the control bandwidth. A constant 5-percent multiplica- 
tive uncertainty is used. The additive uncertainty weight is a high-pass filter which is shaped to envelope 
the frequency response magnitude of each actuator to sensor transfer function above the control band- 
width. To control modes below 1 Hz and gain-stabilize the higher frequency unmodeled or poorly 
modeled modes, the additive uncertainty weight 


W 


add 


0.25 


(* + 5) 2 

(s + 25) 2 


(141) 


is used, as indicated in figure 21. This weight effectively limits the control bandwidth to 1 Hz. The 
diagonal weighting matrix K a(iil provides a separate scaling for the additive uncertainty in each sensor 
channel and is given by K a jj = diag( 1 ,1 ,0.5,25,0.75). The additive uncertainty output associated with 
the MC5 actuator was also scaled by a factor of 0.5. These scalings were determined to provide a tight 
covering for the additive uncertainty weight. The uncertainty structure forD-scaling is 


^ add 0 
® ^ nut It _ 


where is a 5x3 unstructured block and 


'■mult 
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With reference to figure 21, the H x subproblem is defined by 
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(143) 


(144) 


Although several options exist for defining the respective subproblems of the mixed // 2 / H x formula- 
tion, this approach separates performance and robustness according to the most appropriate norms. 

Using this uncertainty structure, D-scales were included in the generalized plant for fixed order 
mixed-norm control design. The generalized plant had 40 states including 26 for the nominal plant, 

6 for the multiplicative uncertainty (a second-order weighting function for each of three plant input 
channels), and 4 each for the left and right frequency varying D-scales. A compensator dimension of five 
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was chosen to result in a true canonical form (the plant has five measured outputs) and reduce 
computation. 

Varying p from 1 to 0.001 in the H 2 subproblem resulted in mixed-norm designs of such low 
authority that the performance was basically the same as open-loop. To allow a higher control band- 
width, the additive uncertainty weight was modified to envelope modes above 3 Hz, but no appreciable 
improvement in control authority was realized. The probable cause of this severely limited control 
authority is the excess conservatism associated with the full unstructured additive uncertainty block. 

Full unstructured uncertainty allows for cross-coupling between inputs and outputs which may not be 
realistic. This uncertainty model also neglects the inherent kinematic coupling of the displacements 
and rates. 

In an attempt to mitigate this excess conservatism, the rate measurements were removed from 
the additive uncertainty, resulting in a 2x3 unstructured additive uncertainty block. However, using addi- 
tive uncertainty weights for both I - and 3-Hz control bandwidths again resulted in similar low authority 
controllers. Thus, with the compensator dimension constrained to five states, unstructured additive un- 
certainty appears to introduce excessive conservatism which severely limits nominal performance in the 
mixed norm design setting. The effect of uncertainty models with compensators having larger dimension 
remains to be assessed. 

Phase II: 

An alternative means of limiting control bandwidth and gain-stabilizing high-frequency modes 
is to use frequency-dependent multiplicative uncertainty in the H ^ subproblem. In this second design 
phase, the additive uncertainty weight is removed and the multiplicative uncertainty at the plant input 
is given by 


^ mult 



(145) 


where W m represents a percentage uncertainty in magnitude of the frequency response at dc. The same 
uncertainty structure for A nm i t given in equation (143) is used for this design phase. 

When using the mixed norm design procedure, it is important to recognize that the mixed norm 
design is only merited when the underlying H 2 controller does not provide robust stability for a given 
uncertainty model and control authority level. Figures 23 and 24 present the robust stability of the full- 
order and reduced-order (5th) baseline Hj designs, respectively, as a function of percent multiplicative 
uncertainty. For less than 10-percent multiplicative uncertainty, the full-order H 2 designs satisfy the 
robust stability constraints and the mixed norm design is not merited. However for 10 percent or larger 
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uncertainty, the full-order H 2 designs are only robustly stable in the low authority range. As expected, 
the shape of the surface in figure 23 indicates that the stability robustness is inversely proportional to 
control authority for a given level of uncertainty and is also inversely proportional to uncertainty at a 
fixed level of control authority. This is a well-known property of H 2 control design that in part moti- 
vates the development of robust control theory. The robust stability plot of the reduced-order H 2 
designs indicates the same trend except that loss of stability due to order reduction precludes high 
authority controllers which violate robust stability to 10-percent uncertainty or larger. The point on 
figure 22 denoted V is the highest authority reduced-order H 2 controller that satisfies robust stability 
for 10-percent multiplicative uncertainty. 



Figure 23. Robust stability of full-order H 2 controllers. 



Figure 24. Robust stability of reduced-order H 2 controllers. 
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Nominal performance of the fixed-order, mixed-norm designs with 5-, I0-. and 20-percent multi- 
plicative (dc) uncertainty is shown in figure 22. Note that for each uncertainty level, a set of mixed-norm 
controllers results with nominal performance in different regions of the cost trade curve. Figure 22 indi- 
cates that the effect of uncertainty is to limit the control authority and, hence, limit the attainable perfor- 
mance. Although this is the same trend observed with the mixed-norm designs that included the additive 
uncertainty model (phase I), the limiting effect of the multiplicative uncertainty is not as severe and con- 
trollers with better nominal performance are obtained. The gaps between the full-order Hi design curve 
and the fixed-order, mixed-norm design curves are indicative of the performance sacrificed to accommo- 
date compensator order constraints. Performance of the mixed-norm designs is severely limited by con- 
straining the compensator to five states. Although the fixed-order, mixed-norm procedure produces 
higher authority stabilizing controllers than does order reduction on the Hi designs, the loss of closed- 
loop stability with Hi order reduction could possibly be remedied with fixed-order // 2 design. 

It should also be noted that although the mixed-norm controllers guarantee robust stability 
according to the uncertainty model for which they are designed, these mixed-norm controllers (except 
for the very low authority designs for 20-percent uncertainty) are not stable with the 40-state evaluation 
model. Thus, the multiplicative uncertainty model does not adequately capture the error introduced to 
the design model during the process of model reduction. Figure 25 presents the frequency response mag- 
nitude of the 26-state design model, the 40-state model, and the envelope of the uncertain design model 
with 10-percent multiplicative uncertainty. Evident from this figure is truncation of modes above 3 Hz 
in the design model as well as the mismatch between models at low frequency that results from model 
reduction. Other input/output frequency responses are similar. 



Figure 25. BLTx to TDSx transfer functions. 
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6. CLOSED-LOOP SYSTEM IDENTIFICATION FOR CONTROLLER REDESIGN 


The traditional approach to control design is to obtain a nominal model of the plant, P , which 
is the basis for control design. Since P is an approximation of the true plant, P, the model based com- 
pensator, Cp, must provide a certain level of robust stability (i.e., Cp must internally stabilize P and 
P ). In the system identification process, an attempt is made to bound the model error P - P , which 
determines the amount of robustness required by the control design. It is also desirable that the control- 
ler provide some level of performance robustness; that is, the achievable performance should not differ 
significantly from the nominal design performance. High-performance control design for a flexible 
space structure is especially challenging since high-fidelity nominal plant models are difficult to obtain. 
The large error bounds that result typically require a very robust, low-performance control design which 
must be tuned on orbit to achieve the required performance. 

An upper bound on achievable performance is 18 > 29 

| ./(/>,C^) | <|| j(P'Cp) 1 + 1 j(p,Cp)-j(p,Cp) | , (146) 

where || j[^P,C p) is the achieved performance, || j{p,Cp ) || is the nominal performance, and 
| ,/^P.C^j - ./(P.Cpj is performance differential. The choice of a specific performance metric ./ and 
norm || * || is determined by the control design methodology. To achieve high performance requires: 

• High nominal performance |<=> || j(^P,Cp j | smallj 

• Robust performance j[p,Cp^J- .l[P,C pj « || j{P,Cpj j . 

The second condition is actually a requirement on model fidelity and indicates that the nominal closed- 
loop system model must closely approximate the actual closed-loop system performance when the com- 
pensator Cp is used. Therefore, the “fitness” of the nominal model is a function of the compensator and 
must be judged from a closed-loop perspective. This fitness is not guaranteed by good open-loop model 
matching nor is it precluded by poor open-loop model matching. 3 

When the model error P - P is large, stability and performance robustness necessitates a low 
authority controller. Since a low authority controller is not as sensitive to model errors, the performance 
differential will be small compared to high authority controllers. However, the performance may not be 
satisfactory. To achieve high performance, the issues of modeling and control design must be treated 
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as a joint problem. The fitness of the design model. P , is a function of Cp. which is itself a function of 
the design plant. In some cases, high-performance control design requires an iterative closed-loop sys- 
tem identification and control design procedure. 

A new closed-loop system identification method is presented in this section which is one step 
of an iterative closed-loop system identification and control design procedure. It is assumed that a mod- 
erately accurate dynamic model of the system to be controlled is available for the initial low authority 
controller design. However, this initial design model is not of sufficient fidelity to permit high authority 
control design. The objective of the closed-loop system identification procedure presented in this chapter 
is to refine the initial control design model based on closed-loop response data. 

In the development of a closed-loop system identification method, consideration must be given 
to the nonuniqueness of the triple (A. in the identified realization. Although there is an infinite 
number of equivalent state space realizations for a system, a system with n states, nu inputs, and ny out- 
puts can be uniquely expressed with a minimum of n(nu + ny) parameters. Having as the objective 
of the closed-loop system identification process the ability to refine an existing design model, one ap- 
proach which circumvents the nonuniqueness problem is to realize the open-loop system matrices 
in a unique, minimal form and directly identify the canonical parameters from closed-loop response 
data. Denery has developed a method of parameter estimation for multivariable state space systems 
from open-loop test data using canonical forms. 43 By utilizing the structure of the closed-loop system 
matrices, an extension of Denery ’s algorithm is developed to estimate the plant parameters based on 
closed-loop response data. 

Proper selection of the objectives of system identification and control design further stresses 
the joint nature of the identification and control problem. Based on the prediction error method, the ob- 
jective of the new closed-loop identification procedure developed in this section is to obtain a model P 
that minimizes the performance differential T — T ^ where T is the actual closed-loop system and T 
is the identified closed-loop system. This system norm cannot be evaluated since T is not known. How- 
ever, the actual and predicted closed-loop measurements are known and an equivalent objective is to 
minimize the prediction error of the closed-loop system, || y — v ||,. The actual and predicted closed-loop 
system outputs, y and y, respectively, are determined for the same set of inputs. The least-squares cost 
functional for control-dependent, closed-loop system identification then is 

J = — [*' (y - y) W{y - y)dt , (147) 

9 J{) 

where IT is a constant matrix chosen to weight the relative importance of different measurement outputs. 
The control objective is matched with the identification objective by designing the controller to mini- 
mize the H~) criterion T 

2 

The rest of the section is structured as follows. First, the algorithm developed by Denery for 
open-loop system identification procedure is presented in detail. Next, the extension of the canonical 
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system identification algorithm to closed-loop system identification is presented. Finally, an algorithm 
for iterative closed-loop system identification and control redesign is presented along with illustrative 
examples. 


6.1 Denerv’s Open-Loop System Identification Algorithm 

In reference 43. a two-step procedure is given which generates parameter estimates based on 
noisy measurement data. The algorithm begins with an equation error procedure, which is similar to a 
linear observer, to generate an initial estimate of the parameters. Noisy measurement data may cause the 
equation error estimates to be biased, but they are sufficiently accurate to initialize the second step, an 
iterative quasi-linearization output error procedure. Since the structure of the two procedures are identi- 
cal with one exception that will be pointed out in the following, these two procedures are combined to 
form an iterative algorithm that is robust to initial parameter estimates and relatively insensitive to 
measurement noise. First, the details of the equation error procedure will be presented, followed by the 
output error procedure. 

6.1.1 Equation Error Procedure 

Consider a model of the state space system to be identified: 


Fz + Git, z(0) = zq 

(148) 

y = Hz . 

(149) 


The objective is to identify some F, G H, and z 0 which represents the dynamics of the unknown 
system based on knowledge of the inputs, u(t), and noisy measurements, y(t). An estimate of the un- 
known parameters may be obtained by minimizing the cost functional given in equation (147). Directly 
minimizing ./ results in a nonlinear optimization process, but in the absence of measurement noise, a 
linear formulation may be obtained. Recognizing that for a perfect model the output in equation (149) 
will exactly equal the measurements and y - v = 0, this difference can be fed back to the model with 
arbitrary gains K and M according to 


: = Fz + G<+ K(y- ht) (150) 

y = Hz+M(y-Hz) , (151) 

which can also be written 

z = F n z + G^u + 8Gi + Ky (152) 

r(0) = -/vo + &o (153) 

v = H n z + My (154) 
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by use of the definitions 


F n = F-KH 

(i 

155) 

G/y = G- 8 G 

(i 

156) 

H n ={1-M)H 

(i 

157) 

1 1 

II 

n 

1 

(] 

'^h 

00 


The elements of (F- KH ) and (I - Mht) may be chosen independently of the unknown elements 
of Fand H by using a maximum of n * ny parameters in K and M when the structure of the system 
is in a specific form. As a consequence. F/y , G j\r, H/y . and "/yo are chosen and the unknown 
parameters are contained in K. M. 8G. and <5 r( ). 

The structure of the system must be such that the unknown parameters in Fare coefficients 
of measured states. To obtain this structure. Denery developed a canonical form for multivariable sys- 
tems which is analogous to a canonical form in reference 62 for multi-input systems. Denery s canonical 
form is called the observer canonical form in reference 63, which is dual to the controller canonical form 
presented in section 3. It can be shown that if the plant dimension is an even multiple of the number 
of outputs and the first n rows of the observability matrix are linearly independent, then the realization 
is canonical and Hwill consist only of ones and zeros. Otherwise, some elements of Hwill be included 
as unknown parameters in the estimation procedure. 

Equation ( 1 54) can now be rewritten as 


y = y N +f{t)y 


(159) 


where y N is the output of the linearized trajectory 


F yZjV + -A'(O) - : N 0 

(160) 

}’N :: H N : N ■ 

(161) 


The sensitivity matrix, fit), is given by 



(162) 
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where the vector of parameters to be estimated is 


vec(K) 

vcc(SG) 

rec(M) 

vec(& {) ) 


( 163 ) 


Using the expression for v from equation (159) in the cost functional, equation ( 147), results in ./ 
becoming quadratic in the unknown parameters. By differentiating./ with respect to the unknown 
parameters and equating the result to zero, the estimate of y is given by 


1-1 r 


I'' f(t) T Wf(t)dt J ( '' f(t) T W{y - y)dt 


or for discrete measurements. 


Y 


im T m) 


S/('f) ! 'M'(.v('/)-v(<,)) 


The /th column of f(t) is Cy, (t), which is the output of the /th sensitivity equation 


(164) 


(165) 
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- F n z y . 


dK 



d{8G) 

dYi 


(166) 



d{Sz q ) 

dYi 


(167) 


y y . - H n z y . + 


dM 

dyd 


(168) 


From y- the system matrices are obtained by solving equations (155H 158) for KG E and r () . These 
values are used in F ^ and as the initial values for the next iteration. 

6.1.2 Output Error Procedure 


If the measurement data used in the equation error procedure are corrupted with unbiased mea- 
surement noise, a bias error will result in the parameter estimates. This can be circumvented by using 
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an output error procedure which yields unbiased parameter estimates based on unbiased noisy measure- 
ments, provided the initial estimates are sufficiently accurate. Typically, the biased estimates obtained 
from the equation error procedure are sufficiently accurate to initialize the output error procedure. 

Hence, the two procedures form a combined algorithm for unbiased parameter estimates based on 
unbiased noisy measurements. 

The output error procedure implements the method of quasi-linearization, which is a well- 
known approach to minimize equation ( 147) subject to equations ( 148M 149). The method of quasi- 
linearization approximates the response of the system model by a nominal trajectory V/y, based on the 
initial parameter estimates, plus a linearized correction about the nominal trajectory. By defining the 
initial estimates F. G H, and r () to be F N , Qy, H N , and r M) , respectively, y may be approximated 
by y A from 

-A - F N~A +[F- Fyyjryy + [Gyy +8G]u, Z A( ) = Z/y () + &() (169) 

>'a = Hn z a + [H ~ Hn] : n • (170) 


where z N is obtained from 


-N = F N z N + GyU. -A'(O) - : N()' >'N ~ ^ N Z N ■ 

Using the definitions in equations ( 155)— ( 156), and recognizing that the initial estimates are sufficiently 
accurate, implies 


KH =K[H n +MH]~KH n (172) 

MH =M[H N +MH]~MH N . (173) 

Substituting these expressions in equations ( 1 69) — ( 1 70) yields 

-A = F n : a +G n u + 8Gi + K y N , z M) -: N q + 8z {) (174) 

y A = H N = A+ My„ , (175) 


which is identical to equations ( 1 52)— ( 1 54) except y N replaces y. Using y A in equation (147) instead 
of y reduces the minimization problem to a form identical to the equation error procedure, except y N 
is used in the place of y when computing/ in the sensitivity equations, equations ( 1 66) — ( 1 68 ). In the 
absence of noise in the measurements, the equation error estimate is the same as the quasi-linearization 
estimate. 
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6.2 Closed-Loop System Identification Algorithm 


Denery's algorithm is extended to closed-loop system identification by expressing the plant in 
observer canonical form and exploiting the structure of the closed-loop matrices. For the plant given by 


v = Av + B\ tv + B 2 u 

(176) 

y = C 2 -V 

(177) 


and a dynamic compensator 


+ 

< 

II 

(178) 

11 = —C c x c , 

(179) 


the resulting closed-loop dynamics are 
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(180) 


"C 2 0 " 


X 

0 I 


_ x c. 


As with the open-loop algorithm, noisy measurements are accounted for by averaging over multiple data 
sets. 


If the plant (A, B 2 -C 2 ) matrices are expressed in observer form, then the C 2 matrix consists 
of ones and zeroes and the unknown parameters in A are coefficients of the measured transformed states. 
(As stated earlier, in some cases the transformation may not be canonical, resulting in the C~> matrix 
having additional nonzero elements, but these parameters can be estimated as well.) It is assumed that 
the compensator state vector time history is recorded. Comparison of equations ( 1 48)— ( 1 49) with 
equations ( 1 80)-( 181) indicates that 


F = 


B c C 2 


-B 2 C c 

A. 


G= 


B\ B 2 
0 0 


and H = 


C 2 0 
0 / 


(182) 


Since the only unknowns in F and G are the plant matrices A, B { , and B 2 , the unknown 
parameter matrices are defined as 


K = 


K 11 
0 


-K v C c 

0 


SG 4 


<5Q| SCj 2 
0 0 


(183) 
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In the case of a true canonical form for the plant, H is completely known. Note that K l j corresponds 
to the unknown parameters in /t, K\j corresponds to the unknown parameters in S 2 - 1 corresponds 

to the unknown parameters in fij, and <5Q 2 = ^\2- 

For closed-loop system identification, the partial derivative expressions in the sensitivity equa- 
tions, equations ( 1 66)— ( 168 ), are modified as follows and the unknown parameter vector is defined as 


rec{K u ) 
vec{K [ 2 ) 
v«r(SQi) 
v('c(& 0 ) 


For Yj corresponding to the (j, k) element of K\ j , 


cJK 

dji 



( 184 ) 


( 185 ) 


and 


— = 0 . — = 0 . ^ = 0 , (18 

dn dYi rjy, 

T 

where eye/ is a matrix of zeros except for a one in the (j, k) element. For y, corresponding to the (j, k) 
element of /C] 2 , 


where 


dK 0 

- 5 — = dy, 

dYi 0 0 


dK 1 2 C c 

dYi 



( 187 ) 


( 188 ) 


which is an nxnc matrix of zeroes except for the jth row which is comprised of the £th row of C c . 
Since <5Q 2 = ^12 > 


d8G 0 

— — = d Yi 
dYi 0 0 


(189) 
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where is zero except for a one in the (j, k) element. The terms and 

are identically zero. For y ,• corresponding to the (j, k) element of £Gj ( , 


dSG 

dy, 


n 

dy, u , 

0 0 


( 190 ) 


where ' ' s zero except for a one in the (j, k) element and all others are identically zero. Similarly, 

r)S~ 

for /, corresponding to the /th element of <5r 0 , is a zero vector with a one in the y'th element 
and all others are identically zero. 


Note that this algorithm is not guaranteed to converge. Since the estimates are determined by 
minimizing the error in the closed-loop time response and not the error in the open-loop plant parameter 
estimates, the plant parameter estimates may not converge to the “true” plant parameters but still provide 
a good control design model. 


6.3 Iterative Control Redesign Examples 

The iterative closed-loop system identification and control design procedure implemented herein 
is patterned after the approach of reference 18 with one notable exception to be pointed out below. 

Iterative closed-loop identification and control redesign algorithm: 

1. Beginning with model P v design a set of // 2 controllers of varying control authority. 

2. Evaluate actual and estimated output and control costs and performance differential 

3. Determine highest performance control design point which satisfies performance differential 

constraint threshold, denoted C ~ .. 

P.i 

4. Using closed-loop data from t|p. Cp . j and P ) , do closed-loop system identification 
to determine P,+\. 

5. Repeat until desired performance is attained. 

This algorithm differs from the framework presented in reference 18 in that the amount of con- 
trol authority increase between iterations is more formally quantified. Recognizing that small changes 
in controller authority tend to result in small changes in performance, a constant scaling factor was used 
in the control design step in reference 18 which was slightly increased each identification/control design 
iteration. Thus the control authority was gradually increased each iteration until an appropriate model 
and high authority control design was achieved. In the procedure introduced above, the performance 
is evaluated for a set of controllers with varying authority to ascertain the onset of performance differen- 
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tial due to model mismatch. Instead of numerous iterations of identification and control design, the 
emphasis is placed on evaluating a set of controllers designed for a common model. By explicitly evalu- 
ating the performance differential for each controller, large steps in control authority may be taken with 
each iteration, resulting in few identification/control design iterations. The output and control costs are 
evaluated in the same manner as was done for the building structure control example of section 5. 
Although the iterative procedure is not guaranteed to converge, the convergence may be checked at each 
iteration by evaluating the performance at each iteration. 

6.3.1 Coupled Mass Example 

As an example of the iterative identification and control procedure, the coupled two mass prob- 
lem illustrated in figure 26 is used. This example problem highlights robust control issues as related to 
flexible space structures and was used as a benchmark problem in reference 64 (with k\ = 0 and d\ - 0 ). 
A disturbance acts on mass two and the control force is applied to first mass. The coefficient matrices are 


0 0 

0 0 

~(k\+k 2 )/ni\ k 2 / niy 
k 2 / m 2 —k 2 / m 2 


0 


0 1 

— (r/| + d 2 ) / ni] d 2 / ni\ 
d 2 / m 2 —d 2 m 2 


(191) 


0 


0 

0 


0 


, b 2 = 

1 / m\ 

0 


1 / m 2 


0 


(192) 


*1 I — *2 



Figure 26. Coupled mass benchmark problem. 
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The slate vector is j - [.v, , .vi , ,vj , .{-> ] T and the measurements are v = [,V| , .vi ] T . so 


10 0 0 
0 10 0 


(193) 


Two cases will be considered— the first having open-loop stable and the second case having a rigid body 
mode. The stable system is described by A] = k 2 = 1.2. m x = m 2 = 1.5, £, = = 0. 1 and the damping 

constant is computed by d = 2C,^/A / m, . 

The procedure begins with an initial plant for control design. As an extreme case, the initial plant 
is obtained by adding 50-percent error to A| , k 2 . £i , and £ 2 and 5-percent error for each of the two 
masses. After transforming the true B 2 ].Q>) triple to observer canonical form, the resulting 

realization is 
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(194) 

0 

0.8000 

0 
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-0.3578 



[51/ 52/] = 
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(195) 


and the corresponding initial triple in observer canonical form is 


A0 = 


0 

0 


-1.1429 

0 

1.1429 ' 


-0.3207 

0 

0.3207 

(196) 

1.1429 

0 

-2.2857 

0.3207 : 

1.0000 

-0.6414 


T 0.6349 
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1 
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(197) 
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In observer canonical form, the C 2 matrix is fixed for a given set of observability indices and the col- 
umns of the A matrix with free elements corresponds to the columns of Cl that have an element equal 
to one. Note that the resulting initial design plant elements varied by 79.28 percent and 42.86 percent 
in the A matrix, and 4.76 percent in the B\ and B 2 matrices from the truth model. 

A set of LQG controllers of varying authority were designed for the initial design plant using 
the weighting matrices 


W = 


0 


0 

Pi nit 


I 


V = 


n 

0 / 


0 


II x 


(198) 


where p is used to vary the control authority. The performance of this set of controllers was then evalu- 
ated with both the design model and the truth model to assess the performance differential that results 
from the initial erroneous model. Recall from the beginning of this section that the performance differ- 
ential is a measure of performance robustness. Figure 27 indicates a large performance difference at all 
control authority levels, so a controller with a moderate authority level (p =5) is chosen for initial 
implementation. 



Figure 27. Performance differential using initial plant: Coupled mass problem. 


Using the LQG controller designed for p = 5 with the initial design model, the closed-loop is 
excited with unit intensity, zero mean random noise low pass filtered at 25 Hz. The closed-loop measure- 
ments are corrupted with a low-intensity, random measurement noise (the standard deviation of the noise 
was equal to 20 percent of the standard deviation of the measurements). In reference 43, measurement 
noise is accounted for by averaging over multiple experiment sets. In this example, five sets of measure- 
ments are used and the five resulting sets of estimated system matrices are averaged. 
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Table 3 gives the initial, actual, and estimated parameters of the system matrices in observer 
canonical form. The significant error is clear as well as the convergence of the parameter estimates after 
50 iterations. As with the combined open-loop algorithm, the first 25 iterations used the equation error 
method and the second 25 iterations used the output error method. The convergence over 50 iterations 
for the correction to the ,4(4,3 ) parameter is shown in figure 28. This parameter corresponds to the 
largest element of y (required the largest correction) at the first correction iteration. A slight discontinu- 
ity is evident at the 25th iteration when the algorithm switched from the equation error method to the 
output error method. However, this is removed after one iteration. 


Table 3. Comparison of initial, actual, and estimated parameters 
for open-loop stable coupled mass problem. 


Initial Parameters 

True Parameters 

Estimated Parameters 

- 1.1429 

- 0.8000 

- 0.8240 

- 0.3207 

- 0.1789 

- 0.1832 

1.1429 

0.8000 

0.8252 

0.3207 

0.1789 

0.1993 

1.1429 

0.8000 

0.8219 

0.3207 

0.1789 

0.1720 

- 2.2857 

- 1.6000 

- 1.6531 

- 0.6414 

- 0.3578 

- 0.3815 

0.6349 

0.6667 

0.6672 

0 

0 

- 0.0064 

0 

0 

- 0.0027 

0 

0 

0.0016 

0 

0 

- 0.0004 

0 

0 

0.0113 

0.6349 

0.6667 

0.6733 

0 

0 

0.0039 



Figure 28. Convergence of AA(4,3). 
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Having refined the initial design model to obtain a more accurate model, a second set of LQG 
controllers is designed and the performance differential evaluated. Figure 29 shows that the gap between 
design performance and achieved performance is considerably decreased at all authority levels when 
compared to figure 27. The identified model results in robust performance (as defined at the beginning 
of this section in regard to performance differential) and good nominal performance. It bears pointing 
out that when the identification experiment was conducted without measurement noise, the achieved 
performance and design performance curves were indistinguishable, indicating that the difference in 
figure 29 is due to measurement noise. More averages and more iterations could possibly further dimin- 
ish the performance differential. Another perspective is to compare the achieved performance for a set 
of LQG designs based on the initial plant model with the achieved performance for the designs based 
on the estimated model, shown in figure 30. The difference between the dashed line and the solid line 
indicates the performance increase attained by performing closed-loop system identification and control- 
ler redesign. Although the amount of performance sacrificed by designing the controller based on an 
initial open-loop design model depends on the specific plant and control design, figure 30 illustrates the 
fact that model error limits achievable performance and better performance can be obtained by refining 
the model to reduce the error. The design point corresponding to p = 5 which was used in the closed- 
loop system identification is indicated by 'x' on figure 30. 



Figure 29. Performance differential using estimated plant: Coupled mass problem. 

A more difficult identification problem is obtained by removing the spring and damper denoted 
k\ and d\ on figure 26, resulting in an unstable rigid body mode in the open-loop plant. Identification 
of open-loop unstable systems (such as spacecraft) is a difficult task which is a further motivation for 
closed-loop system identification. Using the same initial stiffness and damping error, table 4 indicates 
the convergence of the estimated system parameters after 50 iterations. Again, the measurements were 
noise-corrupted and the estimates were averaged after five identification experiments. Similar conver- 
gence of the parameter estimates was observed for this case as for the open-loop stable example. 
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Vu 

Figure 30. Comparison of achieved performance: Coupled mass problem. 


Table 4. Comparison of initial, actual, and estimated parameters 
for open-loop unstable example. 


Initial Parameters 

True Parameters 

Estimated Parameters 

- 1.1429 

- 0.8000 

- 0.8228 

- 0.3207 

- 0.1789 

- 0.1830 

1.1429 

0.8000 

0.8253 

0.3207 

0.1789 

0.1937 

1.1429 

0.8000 

0.8228 

0.3207 

0.1789 

0.1835 

- 1.1429 

- 0.8000 

- 0.8231 

- 0.3207 

- 0.1789 

- 0.1860 

0.6349 

0.6667 

0.6637 

0 

0 

0,0009 

0 

0 

- 0.0011 

0 

0 

- 0.0016 

0 

0 

0.0016 

0 

0 

0.0034 

0.6349 

0.6667 

0.6678 

0 

0 

0.0045 


6.3.2 Building Control Example 

A second example is the building structure used as a control design example in section 5. Using 
the six-state nominal design model as the truth model, the inputs are the ground acceleration and control 
force and the outputs are the relative displacements of the three floors, which results in 30 parameters 
to be estimated. An initial erroneous design model is obtained by adding 5-percent error to the natural 
frequency square terms and the B { and B 2 matrices. Although only 5-percent error is introduced, the 
maximum errors in the elements of the /), By, and B 2 matrices (in observer canonical form) are 

72.2 percent, 9.9 percent, and 39.6 percent, respectively. 
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Using the initial model, a set of LQG controllers is designed using p to vary the control author- 
ity. For this example, filtered noise is used as input excitation and perfect measurements are assumed. 
Figure 31 indicates the performance differential resulting from controllers designed for the initial model, 
which is relatively constant at all control authority levels. A low authority controller is used for closed- 
loop system identification which results in an estimated plant with the performance differential shown 
in figure 32. For this estimated model, there is virtually no performance differential. 



Vu xio -3 

Figure 3 1 . Performance differential using initial plant: Building problem. 



Figure 32. Performance differential using estimated plant: Building problem. 
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6.3.3 Discussion 


Since the end objective is an iterative system identification and model-based control design pro- 
cedure. additional constraints are placed on the input and output processes. The system identification 
is based on closed-loop test data which mandates that the generalized plant, equations (54) — (57), consist 
of actuated inputs and measured outputs only. Figure 33 illustrates this requirement where the disturbance 
inputs, vr and w p , act through the control input and sensor channels and the performance variables. 

" and : p , must be linear combinations of the sensed variables, y, and the control inputs, u. Hence a 
constraint is placed on the generalized plant formulation by the system identification such that the col- 
umns of B { and B p lie in the column space of B 2 and the rows of Cj and C p lie in the row space of C 2 . 
With regard to the H ^ subproblem, this requires the uncertainty representation to consist only of (input 
or output) multiplicative and additive uncertainty models, which implies that it is not possible, for ex- 
ample. to include parametric uncertainty such as mass or stiffness uncertainty in the generalized plant 
formulation. The matrices D 2 \, D 2p , D\ p , and D\ 2 are obtained from the input/output uncertainty 
models and disturbances, and do not contain parameters to be estimated. Consequently, the plant matri- 
ces A. B 2 . and C 2 are the only matrices to be estimated by closed-loop system identification. To relax 
this constraint would require introducing additional actuators and sensors for the sole purpose of system 
identification. 


Disturbances 



Figure 33. Input/output constraint relationships. 

A brief discussion on the closed-loop system identification algorithm from a numerical imple- 
mentation perspective is warranted. The stability of the algorithm is sensitive to several factors. The 
primary factor influencing convergence of the parameter estimation is the matrix inversion in the com- 
putation of y (equation (164) or (165)). For this inverse to exist, the ( ny + nc)xN p sensitivity matrix/ 
must have full row rank where N p = n*(nu + nv + nw) is the number of unknown parameters. Recall 
that the columns of /are time histories of the sensitivity equations, equations ( 166M 168), which include 
the compensator state time histories. If the input is not sufficiently rich to excite the measurement 
and compensator states, then /will not have full rank. This difficulty is compounded as the number 
of parameters increases. If the matrix tends to singularity, the magnitude of the elements of y diverge. 
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In order to alleviate the divergence of y in the examples above, a relaxation factor was introduced that 
scaled y. Scaling y by a relaxation factor of 0.5 typically was sufficient to produce smooth conver- 
gence as seen in figure 28. Without the relaxation factor, the estimates would overshoot and overcorrect, 
resulting in divergence of the parameter estimates. The relaxation factor in essence damped the over- 
shoot of the correction steps at each iteration. This could possibly have been accomplished by using a 
pseudo-inverse of the matrix to zero the small singular values, but that would have introduced error. 
Using the relaxation factor did not introduce error but only slowed the convergence. 
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7. CONCLUSIONS AND RECOMMENDATIONS 


This dissertation has shown that to achieve high performance control requires a control design 
methodology that maximizes performance while guaranteeing robust stability for a given model error. 
The mixed H 2 / H x design procedure presented herein is well suited to this task when used to synthe- 
size a set of controllers that explicitly trade between robustness and performance. Design examples 
have shown that fixed-order, mixed-norm design provides performance comparable to an //-> design 
while overcoming a key deficiency of H 2 design by providing robust stability. When performance 
is defined by an H 2 measure in mixed-norm design, better performance is often obtained than 
the (A -synthesis for the same error bound. Synthesis of the mixed-norm controller set is accomplished 
using a homotopy algorithm which generates a trajectory of gain variations as plant parameters or 
performance specifications vary. 

To improve performance often requires reducing model error through system identification. 

In many cases open-loop testing is not possible, and even when it is, often the most appropriate model 
for control design is obtained from closed-loop response data. A new method for refining a control 
design model from closed-loop response data is presented herein. Based on a prediction error method, 
the open-loop plant parameters are estimated in a canonical form. Examples have shown that higher 
performance can be obtained when the controller is redesigned based on the refined model. 

A major shortcoming of the mixed-norm design procedure is the inherent computational burden 
associated with the homotopy algorithm. Basically a predictor-corrector method, the homotopy algo- 
rithm requires prediction step sizes that are typically quite small due to the ill-conditioning of the Hes- 
sian. In part, the numerical difficulties associated with the often ill-conditioned and indefinite Hessian 
are due to the use of canonical forms for the compensator dynamics. Transforming the controller to 
canonical form results in a more poorly conditioned realization which makes optimization of the control- 
ler parameters more difficult. Future research should investigate the use of over-parameterized realiza- 
tions of the controller, as suggested in reference 49. By using a more robust optimization procedure that 
does not rely on inversion of the Hessian, larger step sizes and more efficient computation could possi- 
bly be attained. The robustness of the continuation algorithm developed herein could possibly be im- 
proved by using newer curve-tracking homotopy algorithms. 65 Another promising approach is to employ 
genetic algorithms which are insensitive to both ill-conditioning of the Hessian and saddle points in the 
search space. To improve the attainable performance, future developments should investigate the incor- 
poration of real parameter uncertainty into the mixed-norm design procedure. Recent advancements in 
mixed real/complex fu synthesis are promising toward this end. 45 Computational efficiency may also 
be gained by taking advantage of developments in computer technology such as parallel architectures. 

Much like the control design homotopy algorithm, a key shortcoming of the closed-loop system 
identification procedure is numerical sensitivity. The solution procedure presented herein requires 
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the inversion of a large data matrix which tends to be ill-conditioned. Ensuring full-rank of the data 
matrix requires sufficiently rich excitation in all of the modes to be identified which is often quite diffi- 
cult in practice due to limitations such as sensor and actuator dynamics. Methods which do not involve 
matrix inversion such as genetic algorithms have potential for the closed-loop parameter estimation 
as well as control design. An additional benefit of genetic algorithms is that the system identification 
can be operating in the background as part of an autonomous identification/controller tuning process. 
Future research should be conducted to that end. Finally, the measurement noise properties (intensity, 
frequency spectrum, etc.) should be explicitly accounted for in the parameter estimation instead of the 
ad hoc use of ensemble averaging. 
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APPENDIX A— NUMERICAL OPTIMIZATION FOR ILL-CONDITIONED 

AND INDEFINITE HESSIANS 


Many numerical optimization problems exist where analytical expressions for the function and 
its first and second derivatives are available. However, numerical determination of the minimum may be 
quite difficult. In particular, ill-conditioned and/or indefinite Hessians present a challenge for numerical 
optimization. A modification to the standard Newton optimization algorithm is presented here that 
effectively accommodates ill-conditioned, indefinite Hessians encountered in the synthesis of fixed-order 
optimal compensators. Synthesis of the optimal control law involves solution of a set of nonlinear 
coupled matrix equations that arise in a minimization problem with explicit first and second derivative 
expressions. However, an accurate initial estimate of the fixed-order compensator gains is often difficult 
to determine due to factors such as order reduction and transformation to certain parameterizations, 
which may result in an ill-conditioned or indefinite Hessian. Before describing the optimization proce- 
dure, preliminaries on quadratic optimization are presented. A more thorough treatment may be found 
in a text on numerical optimization such as references 66 and 67. 

A.l Background 

In many optimization applications, the objective is to determine the parameter vector x 
which minimizes some function /(.v). It is assumed that the function is continuous with analytical expres- 
sions for the continuous first and second derivatives. The gradient vector and Hessian matrix are given 
by 

g(-v) = V/'f a) (199) 


and 


C(.v) = v(v/(.v) r ) = V 2 /(v) . 

respectively, where the gradient operator is 


V(/) 


JL jL , A_ 

d\\ d\2 Sx N 


Necessary conditions for a local minimizes x*, are that 


£(.*-*) = 0 


( 200 ) 


( 201 ) 


( 202 ) 
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and 


s T G(x*)s> 0 se3i N . (203) 

The second order necessary condition, equation (203) is sufficient if 

s T G(.\ *)s > 0 V.v*0 . (204) 

The necessary condition requires the Hessian to be positive semidefinite while the sufficient 
condition requires that the Hessian be positive definite at the local minimum. 

Numerical procedures are used to generate a sequence of points, \ with the objective of 
converging to the fixed solution point .v*. Given an estimate of the solution at the Arth iteration, \ 
a one-dimensional line search in the direction - is used to find the optimum step size, y/^\ which 
minimizes 


/(.v (i+l) ) = /(.v lt) + *< ( *y t )) . (20.‘ 

Typically, the function to be minimized will not be quadratic, but assuming a quadratic approxi- 
mation (or a quadratic ‘model’) of the function is often effective for minimization. From the Taylor 
series expansion of a general function /(.v), 


/(.v + S) = f(x) + g(x) T S + - S t G(.\)8 + o(s T 5) 


(206) 


the quadratic model may be obtained by truncating the terms above second order with the resulting error 


on the order of 
remote from the 


d‘8 


. Near a minimum, quadratic models accurately represent the function. Even when 
minimum, second order information provided by the quadratic model gives insight into 
the most fruitful search directions. Quadratic models also have the property of invariance under linear 
transformation which is important for scaling and preconditioning of the Hessian. 


Based on the Taylor series approximation of a function truncated to second order, Newton's 
method defines the search directions of the quadratic model as 



(207) 


Another means of forming a quadratic model is to utilize search directions that are conjugate 
to the Hessian, given by 


s^ T Gs^= 0 . 


(208) 
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The minimum of a quadratic function can be determined in at most N exact line searches along nonzero 
conjugate directions (Theorem 2.5.1, reference 66). Several options exist to select the conjugate direc- 
tions. but one in particular is discussed below. 

For the minimum of the quadratic model to exist in space, the Hessian must be positive 
definite. However, when is “far" from the minimizer .v*, the Hessian may not be positive definite. 
Ill-conditioning resulting from singularity or near-singularity of the Hessian is a particularly difficult 
situation which must be dealt with by the numerical optimization procedure. The two issues of indefi- 
niteness and ill-conditioning of the Hessian are the motivation for the algorithm presented in the follow- 
ing section. 

A.2 Partitioned Newton Optimization Algorithm 

Ill-conditioning of the Hessian is a key difficulty in minimization of quadratic models. The 
condition number of a matrix is a measure of the closeness of the matrix to singularity and is defined as 
the ratio of the maximum singular value of the matrix to the minimum singular value, where the singular 
values of a matrix G are the square roots of the eigenvalues of G^G . Since the eigenvalues of the Hes- 
sian are a measure of the curvature of the parameter space along directions defined by the corresponding 
eigenvectors, small magnitude eigenvalues indicate very slight curvature, or long “valleys” in the param- 
eter space, which require large steps to achieve the minimum. 

Choosing the eigenvectors of the Hessian as search directions guarantees that the search direc- 
tions are orthonormal and conjugate to the Hessian, thereby forming a quadratic model of the type shown 
in equation (208). 67 When the .v is not “close" to the minimum, the Hessian may be indefinite. In this 
case, Newton's method with one line search in the direction given by equation (207) may not be effec- 
tive. The directions with positive curvature must be searched independent of the directions with negative 
curvature to ensure minimization in each direction of the parameter space. Since the indefiniteness of the 
Hessian indicates that the current estimate, is remote from the minimum, searching along direc- 
tions with negative curvature in the direction of decreasing cost may be effective for moving away from 
critical points such as saddle points. 

An optimization scheme that is robust to ill-conditioned and indefinite Hessians is developed 
which applies Newton’s method to a partitioned Hessian. For the symmetric Hessian G, 

GO = d>A , 

where the eigenvectors, 0y, and eigenvalues. Ay, comprise the matrices 


O - [0i 02 L 0/y] 


(209) 


and 


A — diagjAj, A 2 , L , Ayy} . 


(210) 



The Hessian may be written 


G = OAO r 

= +02^2^2^^ +0/yAyv0A/ (-11) 

N 

= I.<PMl - 
/= 1 


Since the eigenvectors fonn a set of /V conjugate directions, the minimization problem is to 
determine 


where 


T* = arg 


= arg 


nun /|.v^ + i/' Vi +¥^$2 + L +V' (/V Va') 


min f[.x {k) +<t>V 


y = 


|y ; y = 





( 212 ) 


(213) 


This minimization problem consists of N one-dimensional line searches to determine the elements of 

(^*1 . Instead of doing N one-dimensional line searches to determine y *. 


y* = 


If/ O* Xf/ t ' 2 ' 1 '" L 


¥ 


one could approximate the optimal step size in each direction with one scalar, t //, where 


if/ = arg min f 

v 


( N 

x ik) + y/Y<Pi 


V 


/=1 


(214) 


In contrast with the search direction in Newton’s method equation (207), which is a function of the 
inverse of the Hessian, the search direction formed from the Hessian eigenvectors is given by 


^ k) = t<t>i • ( 215 ) 

(=1 

Although insensitive to near-singularity of the Hessian, this approximation of the step sizes would lead 
to similar difficulties as Newton’s method in that directions with large curvature would dominate and 
progress would be restricted to a subspace of the parameter space. 
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A compromise between the approximation of equations (214) and the computational burden 
of the N one-dimensional line searches in equation (212) is to partition the search space into distinct 
subspaces with a line search for each subspace. This approach is implemented in the partitioned Newton 
optimization algorithm, which is formalized below. 

By partitioning the eigenvalue and eigenvector matrices in equations (209)— (2 1 0) into five matrix 


partitions as 


O = [0, o 2 0 3 0 4 q> 5 ] 

(216) 

and 


A = block diagjAj. A 2 ,L .A 3 } , 

(217) 

from equation (21 

1) the inverse of the Hessian may be partitioned as 

4 




(218) 


/=1 


and the Newton search directions analogous to equation (207) may be defined as 

.v (,) = -0, A ~ 1 of g . ( 2|9) 


A.2.1 Algorithm 

1. Given an initial estimate compute the eigenvectors and eigenvalues of the Hessian 
g(.y ( 0) ) and form 0 and A from equations (209 )-(2 10). 

2. Partition 0 and A according to equations (2 1 6) — (2 1 7) where: 

• A] is comprised of “large" positive eigenvalues, 

• A 2 is comprised of “small" positive eigenvalues, 

• A 3 is comprised of “small” negative eigenvalues, 

• A 4 is comprised of “large" negative eigenvalues, 

• A 5 is comprised of “near- zero" eigenvalues, 

and O is compatibly partitioned. The partitions are determined by the intervals 

A/ > A part ^ Ay e -^1 

Ay ol ^ Ay — A part ^ Ay G A 2 
part ^ Ay < ~ ~h(ol Ay G A 3 
Ay < — A p ar t Ay g A4 
— A to l ^ Ay < X lo j => Ay G A3 . 
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To preclude numerical problems due to near zero eigenvalues, A, ro / is defined to effectively set these 
small eigenvalues to zero by neglecting the partition A 5 in the search procedure. 

3. Do four separate one-dimensional line searches to minimize .v according to 

. v ( / ) = .v^'- 1) +^ (, V /) . 

where 

.v (/) =<t>, Aj'o/g, V/ = 1:4 , 

and the sign on is chosen to ensure descent of the cost function. These four searches constitute 
one iteration with ,v^ 4 * used as the initial vector, .r^, for the next iteration. 

Variations on this algorithm may be adopted where line searches are conducted in some of the 
individual directions comprising a partitioned eigenvector matrix, O,, but with increased computation. 
Note that in well-conditioned cases where the Hessian eigenvalues are positive and larger than X part , 
this algorithm reduces to Newton’s method. The following section presents example problems which 
demonstrate the utility of the partitioned Newton algorithm. 

A.3 Examples 

The examples in this section are specifically chosen to highlight ill-conditioned and indefinite 
Hessians. For the line searches, a golden section search is used and Pip ar , = l,A- m , = 10 The stop- 
ping criteria for each test case is based on the relative change in cost between iterations and is given by 



The first two examples are from reference 66 , the third is from reference 67, and the fourth is an optimal 
control design. 

Example 1: 

/(.v) = 1 00( J '2 - vf ) 2 + ( 1 - 1 ) 2 ■ 

This function is known as Rosenbrock’s function and is a well-known test function for optimiza- 
tion methods. Located in a steep-sided parabolic valley, the minimum of this function occurs at 
v * = [i t i] 7 . The Hessian is singular when .v 2 - xf = 0.005, so an extreme test case is to define the 
starting point as V q = [0.0.005]^ - From this initial point, Newton’s method did not improve the cost 
and terminated after one iteration due to the singular Hessian. However, after 13 iterations, the parti- 
tioned Newton algorithm converged to the minimum, accurate to nine decimal places. Figure 34 
graphically depicts the convergence achieved using the partitioned Newton algorithm. 
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1.2 


Minimum 



Figure 34. Convergence of Rosenbrock problem with partitioned Newton optimization. 


Example 2: 

/ (-v) = -vf + -V| A'2 + ( 1 + A‘2 )" . 

For this example, an initial choice is 

• v 0 = [0.0] r => g = [0,2] 7 , G = 

At this starting point, the Hessian is indefinite. The search direction from Newton’s method is 
s = [—2,0] which results in termination since no improvement can be made in the ,v 2 direction. 
The partitioned algorithm breaks this into two orthonormal search directions given by 

A, = 2.4142, .v, = [-0.38268,-0.92388 ] T 

X 2 = -0.4 1 42, s 2 = [0.92388,-0.38268 ] T . 


0 r 

1 2 • 


Table 5 indicates the rapid convergence of the partitioned algorithm to the minimum, even though 
the Hessian is initially indefinite. 


( 220 ) 
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Table 5. Convergence of .v for partitioned Newton method in example 2. 


Iteration 


*2 

0 

0.0 

0.0 

1 

0.6367312144637250 

-1.075920796208131 

2 

0.7011674364048575 

-1 .3467638098297626 

3 

0.6959277688023086 

-1.347963884399868 

4 

0.6958843926891740 

-1.347942196344587 

5 

0.6958843897926642 

-1.347942194896332 


Example 3: 


/(.V) = 100(.v 2 - .vf ) 2 + (I - .V, ) 2 + 90(.v 4 - A ? ) 2 + ( 1 - A , ) : 


+ 10.1 


(_ v 2 — l) +(a*4 — I) -- + 19.8(a‘ 2 — l)(- v 4 — l) • 


( 221 ) 


This example is known as Wood's function and is originally attributed to reference 68. 

Wood’s function is much like Rosenbrock’s function except with four parameters. Beginning at 
A() = [3,-10, -5. 10] r . Newton's method did not converge to the minimum at v* = [l, 1,1,1 ] ! ^ ut terminated 
after 14 iterations at a = [0.83, -2.3 1,-5.32,25.84 ] T ■ Figure 35 indicates that the partitioned algorithm 
successfully found the minimizer a*. The piecewise constant rate of convergence is indicated in figure 
36. It is interesting to note that from v 0 . the Hessian remained significantly indefinite with the minimum 
eigenvalue ranging from — 1.36e— 4 to —87.7 and the partitioned optimization algorithm continued to 
converge to the minimizer. 



Figure 35. Convergence of Wood problem with partitioned Newton optimization. 
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Figure 36. Rate of convergence of Wood problem with partitioned Newton optimization. 


It should be pointed out that with some initial conditions for these functions, Newton's method 
converged more rapidly than the partitioned algorithm, but no case has been observed where the parti- 
tioned method did not converge. 

Example 4: 

The final example demonstrates the partitioned Newton optimization algorithm in the context 
for which it was developed, optimization of fixed-order dynamic compensators. The problem is to design 
an H 2 controller for the coupled mass problem introduced in section 6. Section 3 presents the formula- 
tion of the fixed-order // 2 -optimal control problem with the gradient and Hessian expressions developed 
as a special case of the fixed-order mixed H 2 / H x control problem. This optimization algorithm com- 
prises the local correction step of the homotopy algorithm introduced in section 4. 

For this example, the plant is described by equations (1 91 HI 92) for/) and B 2 , respectively, 
with £| = k 2 = 1 2. wj = m 2 = 1 .5, and £) = £ 2 = 0. 1 . The remaining plant matrices are 

B \ =[/ 0] * Q = [/ of, D u =0, £> 12 -[0 pi] , D 21 = [0 I], D 22 =0 , (222) 


where p is used to vary the control authority. 

In order to illustrate the distinction between the partitioned Newton algorithm and Newton’s 
method, an initial controller is selected to be remote from the minimum. The desired controller 
is tor p = 0.01 and the initial controller is selected as the exact solution for the lower authority case 
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is for p = 0.01 and the initial controller is selected as the exact solution for the lower authority case 
with p = 10. Convergence is defined to occur when the relative change in cost between iterations 
is less than 5xl0 -7 . 

Figure 37 illustrates the convergence of the gradient toward zero with the partitioned Newton 
algorithm as well as the termination of the Newton optimization remote from the minimum. Figure 38 
presents the comparison of the convergence of the (2. 1 ) element of the controller gain matrix. G. to the 
exact value with both the partitioned Newton method and the Newton method. Clearly, the Newton opti- 
mization does not converge to the correct value. For this example, the difficulty in convergence is due 
to the near-singularity of the Hessian which produces large erroneous correction steps with Newton's 
method. These results strongly indicate the deficiencies of an optimization algorithm that uses second 
order information and assumes a positive definite Hessian (such as Newton’s method). The partitioned 
Newton algorithm presented herein is better suited for these problems by accounting for ill-conditioned 
and indefinite Hessians. 

Norm of Gradient 



Figure 37. Convergence of gradient for control design. 
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Figure 38. Convergence of control gain estimate. 


The initial and exact final controller gains for the controller in the controller canonical form 
(defined in equation (35)) are 


0.4073 

-0.6215 

-0.0003 -0.3032' 



Gn = 

2.3727 

4.1798 

4.3888 2.2904 


(223) 

T 8.3899 

-37.0769 

-1.2240 -13.8905' 



exact — 29.5261 

40.7987 

27.6252 9.8040 

• 

(224) 

Using the partitioned Newton method, after 

1 1 iterations the percent error in each element is 



G error 


0.0562 

0.0637 


0.0650 

0.0613 


0.1305 

0.0608 


0.0671 

0.0637 


percent . 


The Newton method terminated after 18 iterations with percent error 


G crro 


224 

215 

698 

195 

223 

253 

245 

275 


percent . 


(225) 


(226) 
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